NPS71-88-001    (Second  Revision) 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


A  PROGRAM  TO  COMPUTE  MAGNETIC  ANOMALY 
DETECTION  PROBABILITIES 


R.  N.  FORREST 

October  1988 

(Revised  March  1990) 


FedDocs 
D   208.14/2 
NPS-7 1-88-00  V990 


Approved  for  public  release;  distribution  is  unlimited. 

Prepared  for: 

Naval  Postgraduate  School 
iterey,  CA  93943-5000 


i^r^w-^-^0  cx 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


Rear  Admiral  R.  W.  West,  Jr. 
Superintendent 


Harrison  Shull 
Provost 


This  report  was  prepared  in  conjunction  with  »-^"»* UCted 
for  the  Chief  of  Naval  Operations  and  funded  by  the  Naval 
Postgraduate  School: 


UNCLASSIFIED 


ECURITY  CLASSIFICATION  OF  THIS  PAGE 


°^LEY  KNOX  LIBRARY 
NAVAL  POSTGRADUATE  SCHno. 
MONTEREY  r,  r^fC 


m- 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


a    REPORT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


lb    RESTRICTIVE  MARKINGS 


!a.  SECURITY  CLASSIFICATION  AUTHORITY 


!b.  DECLASSIFICATION /DOWNGRADING  SCHEDULE 


3     DISTRIBUTION /AVAILABILITY  OF  REPORT 

Approved    for   public    release;    distribution 
is    unlimited. 


I    PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 


NPS71-88-001    (Second    Revision) 


5.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 


)a.  NAME  OF  PERFORMING  ORGANIZATION 

Naval    Postgraduate    School 


6b    OFFICE  SYMBOL 
(If  applicable) 

AW 


7a.  NAME  OF  MONITORING  ORGANIZATION 

Naval    Postgraduate    School 


ic  ADDRESS  (City.  State,  and  ZIP  Code) 

Monterey,    CA      93943-5000 


7b    ADDRESS  (City.  State,  and  ZIP  Code) 

Monterey,    CA      9394  3-5000 


la.  NAME  OF  FUNDING  /SPONSORING 
ORGANIZATION 

Naval  Postqraduate  School 


8b    OFFICE  SYMBOL 
(If  applicable) 


9    PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

Direct    Funding,    O&MN 


Ic  ADDRESS  (City,  State,  and  ZIP  Code) 


Monterey,  CA.  93943 


10    SOURCE  OF  FUNDING  NUMBERS 


PROGRAM 
ELEMENT  NO 


PROJECT 
NO 


TASK 
NO 


WORK  UNIT 
ACCESSION  NO 


1.  TITLE  (Include  Security  Classification) 

A  PROGRAM  TO  COMPUTE  MAGNETIC  ANOMALY  DETECTION  PROBABILITIES 


2.  PERSONAL  AUTHOR(S) 


R.  N.  FORREST 


3a.  TYPE  OF  REPORT 

Technical 


13b    TIME  COVERED 
FROM  TO 


14.  DATE  OF  REPORT   (Year,  Month.  Day) 

October      1988 


15.  PAGE  COUNT 

70 


6.  SUPPLEMENTARY  NOTATION 


Revised  Mar  1990) 


COSATi  CODES 


FIELD 


GROUP 


SUB-GROUP 


18   SUBJECT  TERMS  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

Magnetic  Anomaly  Detection,  MAD,  Lateral  Range  Function, 
Detection  Models,  Magnetic  Dipole 


1 9  ABSTRACT  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

The  report  contains  user  instructions,  a  listing  and  documentation  for  a  microcomputer 
BASIC  program  that  can  be  used  to  compute  an  estimate  of  the  probability  that  a  magnetic 
anamoly  detection  (MAD)  system  such  as  the  AN/ASQ-81  will  detect  a  submarine  during 
an  encounter. 


20.  DISTRIBUTION /AVAILABILITY  OF  ABSTRACT 

□  UNCLASSIFIED/UNLIMITED       □  SAME  AS  RPT  □  DTIC  USERS 


21 


ABSTRACT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


22a    NAME  OF  RESPONSIBLE  INDIVIDUAL 
R.    N.     FORREST 


22b  TELEPHONE  (Include  Area  Code) 

(408)    646-2653 


22c.  OFFICE  SYMBOL 
OR/Fo 


>D  Form  1473,   JUN  86 


Previous  editions  are  obsolete 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


The  second  revision  of  the  report  contains  corrections  to  some 
errors  that  were  contained  in  the  original  version  as  well  as  some 
that  were  introduced  in  the  first  revision.  In  addition,  on  Page 
43  there  is  a  discussion  of  the  definition  of  sweep  width  as  it 
relates  to  the  lateral  range  functions  that  are  generated  by  the 
detection  models  on  which  the  program  is  based.  As  was  noted  in 
the  earlier  versions  of  the  report,  the  program  in  the  report  is 
presented  without  representation  or  warranty  of  any  kind. 


Table  of  Contents 

I .  Introduction 1 

II .  Program  User  Instructions 3 

III .  Encounter  Model  Limitations 6 

IV.  Parameter  Values 9 

Appendix  1.   The  Detection  Statistics 12 

Appendix  2.   Program  Probability  Calculations 18 

Appendix  3 .   The  Magnetic  Signal 20 

Appendix  4 .   The  Anderson  Formulation 24 

Appendix  5.   The  Encounter  Equations  of  Motion 3  0 

Appendix  6.   The  Submarine  Magnetic  Dipole 3  3 

Appendix  7.   The  Earth  Magnetic  Field  and  Dip  Angle • •  3  6 

Appendix  8.   An  Alternative  Encounter  Model 41 

Appendix  9.   An  Example  of  the  Program  Output 4  3 

Appendix  10.  A  Program  Listing 54 

References 61 


I.    Introduction 

This  report  contains  user  instructions,  a  listing  and 
documentation  for  a  microcomputer  BASIC  program  that  can  be  used 
to  compute  an  estimate  of  the  probability  that  a  magnetometer 
based  detection  system  such  as  the  AN/ASQ-81  will  detect  a 
submarine  magnetic  dipole  field  during  an  encounter. 

The  program  generates  detection  probabilities  based  on  two 
encounter  models.   In  the  first  encounter  model,  the  detection 
system  uses  a  crosscorrelation  detector.   In  the  second  encounter 
model,  the  detection  system  uses  a  square  law  detector.   Relative 
to  operationally  realizable  values,  probabilities  based  on  the 
first  model  represent  upper  bounds  and  those  based  on  the  second 
model  represent  lower  bounds.   For  both  encounter  models,  the 
signal  is  proportional  to  the  square  of  the  magnitude  of  the 
projection  of  the  dipole  field  on  the  earth  magnetic  field  and 
the  magnetic  noise  does  not  change  with  changes  in  the  position 
of  the  magnetometer.   Also,  for  both  encounter  models,  an 
encounter  is  a  straight  line  encounter  with  constant  vertical 
separation. 

The  encounter  models  can  be  interpreted  as  models  of  a 
magnetic  anomaly  detection  (MAD)  system  on  an  aircraft  that  is 
moving  with  constant  course,  speed  and  altitude  in  an  encounter 
with  a  submarine  moving  with  constant  course,  speed  and  depth. 
Or,  they  can  be  interpreted  as  models  of  a  stationary  magnetic 
anomaly  detection  system  in  an  encounter  with  a  submarine  moving 
with  constant  course,  speed  and  depth. 


The  program  parameters  include:   encounter  latitude  and 
longitude,  submarine  induced  magnetic  moments,  submarine 
permanent  magnetic  moments,  submarine  course,  speed  and  altitude, 
magnetometer  course,  speed  and  depth,  encounter  lateral  range 
(the  horizontal  range  at  the  closest  point  of  approach  in  a 
straight  line  encounter)  and  false  alarm  rate.   (In  the  program, 
a  false  alarm  is  the  event  that  the  detection  system  classifies 
noise  as  a  dipole  signal.) 


II.   Program  User  Instructions 

As  listed  in  Appendix  10,  the  program  can  be  run  under  a 
BASIC  language  that  is  compatible  with  IBM  PC  BASIC.   If  the 
listing  is  used  to  enter  the  program  through  the  keyboard,  then 
the  program  should  be  saved  with  the  name  MAD.BAS  or  the  value  of 
N$  on  line  40  should  be  changed  to  the  file  name  under  which  it 
is  saved. 

The  program  contains  user  instructions  in  the  form  of  query 
and  parameter  limitation  messages.   As  an  example  of  the  former, 
after  starting  the  program  under  BASIC,  the  following  message 
should  appear: 

Magnetic  Anomaly  Detection  (MAD)  Lateral  Range  Function 

generate  or  print  a  program  data  file  (g/p) ? 

By  entering  p,   data  can  be  printed  from  a  program  data  file 
that  was  generated  by  the  program.   By  entering   g,  a  program 
data  file  can  be  generated  for  a  set  of  user  specified 
conditions.   With  either  response,  a  sequence  of  additional 
queries  is  displayed.   These  queries  require  either  an  indicated 
response  or  a  parameter  value  as  the  input.   If  the  initial 
response  is   g,  the  sequence  includes  queries  whose  responses 
determine  whether  or  not  an  auxiliary  data  file  will  be  generated 
that  can  be  used  for  future  input  of  magnetic,  processing  or 
kinematic  data.   In  particular,  the  first  query  in  the  sequence 
gives  the  option  of  using  a  combined  magnetic,  processing  and 


kinematic  data  file.   If  the  response  indicates  that  it  should  be 
used  and  the  file  is  available,  the  parameter  values  that  remain 
to  be  entered  in  order  to  generate  a  program  data  file  are  the 
following:  the  false  alarm  rate,  the  magnetic  noise,  the  maximum 
encounter  lateral  range  and  the  lateral  range  step.   The  combined 
file  should  be  used  only  if  the  effect  of  varying  just  one  or 
more  of  these  parameters  is  desired.   If  the  response  indicates 
that  the  file  will  not  be  used,  queries  concerning  magnetic, 
kinematic  and  processing  parameter  values  are  displayed. 

After  all  of  the  program  parameter  values  have  been  entered, 
a  query  is  displayed  giving  the  option  of  generating  the  combined 
file.   Then  a  query  giving  the  option  of  generating  a  program 
data  file,  a  query  giving  the  option  of  printing  the  encounter 
parameter  values  and  a  query  giving  the  option  of  printing 
lateral  range  function  values  are  displayed.   The  lateral  range 
function  values  are  the  encounter  detection  probabilities  indexed 
by  lateral  range.   The  parameters  maximum  lateral  range  and 
lateral  range  step  determine  the  index  lateral  ranges  of  the 
encounters  for  which  probabilities  of  detection  are  computed. 

The  program  generates  magnetic  signal  values  that  correspond 
to  points  in  time  during  an  encounter.   Following  the  lateral 
range  function  query,  a  query  is  displayed  that  gives  the  option 
of  printing  the  magnetic  signal  values  for  an  encounter.   If  the 
option  is  exercised,  the  option  is  repeated.   When  the  option  is 
not  exercised,  a  query  is  displayed  that  gives  the  option  of 


generating  or  printing  a  new  program  data  file.   If  this  option 
is  not  exercised,  the  program  ends. 

Some  suggested  guides  for  determining  parameter  values  can 
be  found  in  Section  IV  of  this  report. 


III.  Encounter  Model  Limitations 

In  the  two  encounter  models,  the  signal  is  the  unfiltered 
signal  that  is  generated  by  a  magnetic  dipole  that  moves  relative 
to  a  magnetometer.   Depending  on  the  input  filtering,  describing 
a  ship  magnetic  anomaly  as  a  dipole  anomaly  should  not  be  a 
significant  limitation  for  an  encounter  slant  range  at  the 
closest  point  of  approach  (CPA)  that  is  greater  than  one  hull 
length.   The  detection  decision  is  based  on  samples  from  a  single 
time  interval  (window)  that  is  centered  on  the  CPA.   The  length 
of  the  sampling  interval  and  the  sampling  rate  are  parameter 
values  that  are  inputs  to  the  program.   In  terms  of  signal-to- 
noise  ratio,  there  is  an  optimum  sampling  interval  length 
(integration  time)  and  sample  rate.   Although  dipole  signal 
energy  is  not  symmetrically  distributed  about  the  CPA  time,  for  a 
given  sampling  interval,  the  difference  between  the  signal  energy 
for  the  optimum  interval  location  and  the  CPA  centered  location 
should  not  be  significant  in  most  cases. 

The  encounter  model  magnetometer  noise  samples  are 
determined  by  a  gaussian  random  process.   They  are  values  of 
independent  identically  distributed  normal  random  variables.   The 
standard  deviation  of  these  random  variables  is  referred  to  as 
the  magnetic  noise  and  the  variance  is  the  noise  in  the  sense  of 
the  signal-to-noise  ratio.   The  degree  of  correspondence  between 
this  process  and  operational  noise  depends  on  the  nature  of  the 
dominant  operational  noise  sources  and  on  the  magnetometer  input 
filtering  (noise  whitening) . 


In  the  encounter  models,  the  intervals  are  adjacent  but  not 
overlapping  and  a  detection  statistic  corresponds  to  each  sample 
interval  and  its  value  is  determined  by  the  sample  values.   If 
the  value  of  the  detection  statistic  for  an  interval  equals  or 
exceeds  a  threshold  value,  a  detection  is  indicated,   The 
threshold  value  is  determined  by  the  false  alarm  probability 
which  in  turn  is  determined  by  the  false  alarm  rate  and  the 
sampling  interval  length. 

A  detector  that  used  a  moving  sample  interval  that  was 
generated  by  replacing  the  oldest  sample  by  the  newest  one  would 
correspond  more  closely  to  the  detector  in  an  operational 
detection  system.   In  a  model  of  the  detector,  since  the  sample 
windows  overlap,  the  detection  statistics  would  represent  a 
sequence  of  random  variables  that  were  correlated  over  an 
interval  equal  to  the  width  of  the  sample  interval.   Because  of 
this  dependence,  it  seems  unlikely  that  the  results  that  would  be 
obtained  with  an  encounter  model  based  on  an  overlapping  interval 
detector  would  differ  significantly  from  those  obtained  with  an 
encounter  model  based  on  a  nonoverlapping  interval  detector. 

For  encounter  lengths  of  the  order  of  a  few  nautical  miles 
or  less,  the  straight  line  encounter  condition  should  not  be  a 
significant  limitation.   In  particular,  this  should  be  the  case 
for  a  fixed  magnetometer  since,  for  a  submarine  (or  surface  ship) 
target,  vertical  separation  and  course  changes  should  be  less 
likely  to  occur. 


Other  models  are  available  that  can  be  used  as  the  basis  for 
computing  an  estimate  of  the  probability  that  a  magnetic  anomaly 
detection  system  will  detect  a  submarine  during  an  encounter. 
For  example,  one  that  is  described  in  Appendix  8  can  be  used  to 
determine  the  slant  ranges  of  straight  line  encounters  for  which 
the  detection  probability  is  equal  to  a  specified  value.   The 
parameter  values  that  are  required  to  do  this  are  an  average 
submarine  dipole  moment,  a  detection  system  capability  factor  and 
a  noise  factor.   Values  for  these  parameters  can  be  determined 
from  operational  data.   However,  the  values  are  specific  to 
averages  over  a  particular  set  of  encounter  conditions.   An 
advantage  of  the  two  encounter  models  relative  to  this  model  is 
their  adaptability  to  different  magnetic,  processing  and 
kinematic  conditions. 


IV.   Parameter  Values 

The  magnetic  parameter  value  queries  are  generally 
explanatory  with  regard  to  the  value  that  should  be  entered. 
This  is  also  true  of  the  kinematic  parameter  value  queries. 
However,  there  is  some  ambiguity  with  respect  the  processing 
parameter  value  queries  and  the  noise  parameter  value  query.   To 
reduce  this  ambiguity,  a  brief  discussion  of  the  common 
characteristics  of  the  two  encounter  model  processing  parameters 
and  the  noise  parameter  is  given  below.   This  is  followed  by  some 
guidelines  for  choosing  these  parameter  values. 

In  both  encounter  models,  a  decision  is  made  at  the  end  of 
each  sampling  interval.   The  decision  is  either  noise  energy  was 
present  during  the  interval  or  noise  energy  plus  signal  energy 
was  present  during  the  interval.   The  sample  intervals  are 
adjacent,  equal  width,  nonoverlapping  time  intervals.   The  number 
of  samples  that  are  input  in  a  sample  interval  is  determined  by 
the  sampling  rate  and  the  interval  length. 

The  program  default  choice  for  the  sampling  rate  is  2 -MAXF 
where  MAXF  is  a  parameter  that  is  labeled  the  maximum  magnetic 
signal  frequency.   This  sample  rate  is  the  Nyquist  rate  for  an 
ideal  low  pass  filter.   However,  the  signal  in  a  sample  interval 
that  is  computed  by  the  program  represents  an  unfiltered  dipole 
signal.   This  is  a  reasonable  approximation  if  the  signal  energy 
that  is  associated  with  signal  components  greater  than  MAXF  is 
relatively  small.   As  discussed  below,  the  noise  energy  in  a 
sample  interval  should  be  considered  to  be  proportional  to  MAXF 


in  order  to  be  consistent  with  the  encounter  models.   Ideally,  a 
default  choice  for  MAXF  would  make  the  ratio  of  the  signal  energy 
to  the  noise  energy  a  maximum  for  the  sampling  interval  of  an 
encounter.   The  program  default  choice  for  MAXF  is  2 -MAXVM/MINRO 
where  MAXVM  is  a  user  estimated  maximum  encounter  relative  speed 
converted  from  knots  to  meters  per  second  and  MINRO  is  a  user 
estimated  maximum  slant  range  at  CPA  in  meters  in  terms  of  a  just 
detectable  target.   (A  precise  definition  of  a  just  detectable 
target  can  be  made  in  terms  of  a  specified  detection  probability, 
false  alarm  probability  and  target  dipole  moment.)   The  default 
choice  for  MAXF  is  consistent  with  the  observation  in  Reference  1 
that  if  an  optimum  value  for  MAXF  is  determined  for  a  minimum 
dipole  moment  target,  then  no  significant  increase  in  MAXF  is 
required  in  order  to  maintain  a  required  detection  probability  if 
the  encounter  lateral  range  is  decreased  even  though  the  signal 
energy  spectrum  is  shifted  to  higher  frequencies. 

Because  of  the  detection  statistics  and  the  gaussian  noise 
model,  if  there  were  no  penalty  for  decision  delay,  a  sample 
interval  length  for  a  signal  should  be  chosen  equal  to  the  signal 
duration,  since  this  would  make  the  detection  probability  for  an 
encounter  a  maximum.   In  the  program,  a  default  choice  for  the 
sample  interval  length  (integration  time)  is  2 -MAXRO/MINVM  where 
MINVM  is  a  user  estimated  minimum  encounter  relative  speed 
converted  from  knots  to  meters  per  second  and  MAXRO  is  a  user 
estimated  maximum  slant  range  at  CPA  in  meters  for  a  detectable 
encounter  (in  terms  of  a  specific  detection  probability  and  false 
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alarm  probability)  with  a  user  specified  maximum  dipole  moment 
target.   This  choice  might  be  considered  a  balance  between 
minimizing  decision  delay  and  maximizing  detection  probability. 
In  the  encounter  models,  the  sample  intervals  are  located  so  that 
the  CPA  time  is  at  the  center  of  a  sample  interval  and  this 
sample  interval  is  the  only  one  that  contains  signal  energy. 
This  characteristic  is  consistent  with  the  default  choice  for  the 
sample  interval.   For  a  given  sample  interval  length,  although  in 
general  the  interval  location  is  not  the  optimum  one  in  terms  of 
signal  energy,  it  should  be  approximately  so  in  most  cases. 
The  program  noise  parameter  is  SIG.   It  represents  the 
standard  deviation  a      associated  with  the  magnetic  noise  process 
of  the  two  encounter  models.   In  terms  of  the  ideal  low  pass 
filter  implied  by  the  encounter  models,  its  square  should  be 
equal  to  MAXF-(SIGO)2  where  SIGO  is  the  magnetic  noise  process 
constant  power  spectral  density.   The  program  does  not  enforce 
this  relation.   Therefore,  in  using  the  program,  the  implied 
relation  between  the  two  input  parameters:  magnetic  noise  and 
maximum  magnetic  signal  frequency  should  be  kept  in  mind.   If  an 
average  value  of  the  peak-to-peak  magnetic  noise  for  an  encounter 
can  be  estimated,  for  example  from  a  magnetometer  trace,  then  the 
value  for  SIG  should  be  chosen  so  that  the  estimate  is  4  to  6 
times  this  value. 
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Appendix  1.   The  Detection  Statistics 

In  this  appendix,  y,,  y2,  •  •  •  ,  yM  are  sequential  values 
(voltages)  representing  the  sample  values  in  a  sample  interval. 
They  are  the  input  to  a  magnetometer's  detector.   With  these  M 
sample  values,  the  detector  computes  the  value  of  a  detection 
statistic.   This  value  is  represented  by  x  and  the  detection 
decision  corresponding  to  the  sample  interval  is  determined  by 
the  decision  rule:  If  x  >  x*,  then  the  input  during  the  sample 
interval  was  noise  plus  signal,  otherwise,  the  input  was  noise. 

For  both  encounter  models,  the  detection  probability  and  the 
false  alarm  probability  are  decreasing  functions  of   x*   and  the 
relation  is  one-to-one  in  both  cases.   In  the  program,  the  false 
alarm  probability  pF   is  used  to  determine  a  unique  value  of 
the  threshold  x*.   This  value  is  then  used  to  determine  a  unique 
value  of  the  detection  probability  pD. 

In  the  program,   pF   is  found  using  the  relation  pF  =  R-6t 
where  R   is  the  false  alarm  rate  in  false  alarms  per  second 
and  6t      is  the  sample  interval  length  in  seconds.   This  relation 
is  based  on  the  following  argument:   With  no  signal  energy  in  a 
sample  interval,   y,  =  n, ,  y2  =  n2,  •  •  •  ,  y  =  nM  where  n,,  n2, 
,  nM  are  noise  values  (voltages)  input  to  a  magnetometer's 
detector  that  are  values  of  independent  normal  (gaussian)  random 
variables,  each  with  mean  zero  and  standard  deviation  a. 
Because  of  this,  in  the  encounter  models,  values  of  x   for 
different  sample  intervals  are  the  values  of  independent  random 
variables  that  determine  two  outcomes:   x  >  x*  or  x  <  x*. 
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Therefore,   pF   is  the  same  for  each  sample  interval  and,  in 
terms  of  these  outcomes,  the  sequence  is  a  series  of  repeated 
independent  Bernoulli  trials  and  the  expected  number  of  trials 
between  false  alarms  is   l/pF.   Since  the  time  between  trials  is 
«5t,  the  expected  number  of  seconds  between  false  alarms  is   <St/pF 
and  the  expected  number  of  false  alarms  per  second  is  equal  to 
pF/6t.   So,  R  =  pF/<5t. 

The  determination  of  x*  depends  on  the  encounter  model 
statistic.   For  both  encounter  models,  when  there  is  a  signal, 
y,  =  n,  +  s,,  y2  =  n2  +  y2,  •  •  •  ,  yM  =  nM  +  sM  where   s,,  s2,  •  •  ■  , 
sM   are  signal  values  (voltages)  input  to  a  magnetometer's 
detector.   The  models  imply  that  the  signal  values   s  =  c-Hs 
where  c   is  a  constant  whose  value  is  determined  by  the 
characteristics  of  the  encounter  magnetometer  and  where  the  Hs 
are  dipole  magnetic  signal  intensities.   The  models  also  imply 
that  the  noise  values   n   are  determined  by  a  gaussian  stochastic 
processes  characterized  by  a  standard  deviation  a      and  that 
n  =  c-HN  +  n1   where  the  HN  are  magnetic  noise  intensities  and 
the   n'   are  magnetometer  instrument  noise  values  (voltages).   In 
the  program,  the  magnitude  of   c   is   1.   Since,  for  both  models, 
pD  depends  only  on  the  ratio  of  signal  energy  to  noise  energy 
for  a  sample  interval,  this  is  a  satisfactory  choice  for  the 
program  if  instrument  noise  process  is  assumed  to  be  independent 
of  the  magnetic  noise  process  and  to  be  expressed  in  terms  of  an 
equivalent  magnetic  noise  H  =  (1/c)  n'. 
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For  both  encounter  models,  the  signal  (the  average  signal 
power)   S  =  (1/M) -2  sK2  and  the  noise  (the  expected  value  of  the 
average  noise  power)   N  =  a2       so  that  the  signal-to-noise  ratio 
is   (1/M)  -S  sK2/cr2  where  the  sum  index   K  =  1,  2,  •••  ,M. 
The  Crosscorrelation  Detector  Statistic:   The  statistic  for  the 
first  encounter  model  is  a  crosscorrelation  detector  statistic 
that  is  defined  by  the  sum 

x  =  S  yK-sK 
where  the  summation  index  K  =  1,  2,  • • •    ,  M  and  the  sum  is  over 
the  values  corresponding  to  a  sample  interval.   For  the  first 
encounter  model,  the  characteristics  of  both  the  noise  and  the 
signal  are  reguired  in  order  to  determine  encounter  detection 
probabilities.   In  particular,  the  signal  values  for  an  encounter 
are  in  the  memory  of  the  detector  prior  to  the  encounter.   For 
the  encounter  conditions  and  a  specified  false  alarm  probability, 
the  statistic  is  optimum  in  the  sense  that  the  encounter 
detection  probability  for  this  statistic  is  at  least  egual  to 
that  for  any  other  statistic.   Because  of  these  considerations, 
encounter  probabilities  based  on  the  crosscorrelation  statistic 
can  be  considered  to  represent  upper  bounds  on  detection 
performance  against  dipole  targets  for  magnetometers  of  the  type 
described  by  the  models. 

For  a  sample  interval  without  signal  energy,  x   is  the  value 
of  a  normal  random  variable  with  a  mean  nx  =   0   and  a  variance 
ax2  =  a2 -2  sK2  where  a   is  the  standard  deviation  associated  with 
the  noise  process  and  the  sum  index   K  =  1,  2,  • • •  ,  M  and  the 
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sum  is  over  the  values  corresponding  to  the  sample  interval. 
This  implies  that 

pF  =  1  -  P(x'/ax) 
where   P(z)   is  the  standard  normal  cumulative  distribution 
function.   This  relation  is  the  basis  for  determining  the 
threshold  value   x*. 

For  the  sample  interval  with  signal  energy,  x   is  the  value 
of  a  normal  random  variable  with  a  mean  /xx  =  2  sK2  where  the  sum 
index   K  =  1,  2,  • • •  ,  M  and  the  sum  is  over  the  values 
corresponding  to  the  sample  interval.   This  implies  that 

pD  =  1  -  P(v*  -  d1/2) 
where  v*  =  x'/ax   and   d  =  2  s*/o2  =     M-S/N.   This  relation  is  the 
basis  for  determining  encounter  detection  probabilities  for  the 
first  encounter  model.   The  relation  implies  that  for  a  specified 
false  alarm  probability   pF   the  detection  probability   pD  is  an 
increasing  function  of  the  signal  to-noise  ratio   S/N. 
The  Square  Law  Detector  Statistic:   The  statistic  for  the  second 
encounter  model  is  a  square  law  (energy)  detector  statistic  that 
is  defined  by  the  sum 

x  =  E  yK2 
where  the  sum  index   K  =  1,  2,  • • •  ,  M   and  the  sum  is  over  the 
values  corresponding  to  the  sample  interval.   For  the  second 
encounter  model,  only  the  characteristics  of  the  noise  are 
required  to  determine  encounter  detection  probabilities. 
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For  a  sample  interval  without  signal  energy,  x/a2   is  the 
value  of  a  chi-sguare  random  variable  with  M  degrees  of 
freedom.   This  implies  that 

pF  -  1  -  P(x*/cj2|M) 
where   P(x'/a2|M)   is  the  chi-sguare  cumulative  distribution 
function  for  a  chi-sguare  random  variable  with  M  degrees  of 
freedom  and  where  a      is  the  standard  deviation  associated  with 
the  noise  process.   This  relation  is  the  basis  for  determining 
the  threshold  value   x*. 

For  the  sample  interval  with  signal  energy,  x/a2   is  the 
value  of  a  noncentral  chi-sguare  random  variable  with  M  degrees 
of  freedom  and  noncentral  parameter  2  sK2/a2  where  the  sum  index 
K  =  1,  2,  • • •  ,  M  and  the  sum  is  over  the  values  corresponding 
to  the  sample  interval.   This  implies  that 

pD  =  1  -  P(x*/a2|M,  2  sK2/a2) 
where   P(x'/a2|M,  E  sK2/a2)   is  the  cumulative  distribution 
function  for  a  noncentral  chi-sguare  random  variable  with  M 
degrees  of  freedom  and  noncentral  parameter   2  sK2/a2=  M-S/N. 
This  relation  is  the  basis  for  determining  encounter  detection 
probabilities  for  the  second  encounter  model.   The  relation 
implies  that  for  a  specified  false  alarm  probability  pF,  the 
detection  probability  pD  is  an  increasing  function  of  the 
signal-to-noise  ratio  S/N.   This  is  made  more  evident  by  the 
following  relation: 

P(x*/a2|M,  2  sK2/a2)  =  2{  (aJ/Ji )  -exp(-a)  -P[x*/a2|  (M  +  2  •  J)  ]  } 
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where  the  parameter   a  =  (1/2)  -2  s*/o2   =  (M/2)  •  (S/N)  ,  the  sum 

index  J  =  0,  1,  2,  • • •   and  the  sum  index   K  =  1,  2,  • • ■  ,  M, 

Note,  since   P(x*/a2|M)  <  P[x'/a2|(M  +  2  •  J)  ]   for  J  =  0,  1,  2, 

,  the  relation  pD  >  pF   is  satisfied. 
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Appendix  2.   Program  Probability  Calculations 

The  program  evaluates  the  cumulative  and  inverse  cumulative 
distribution  functions  using  approximations  described  in 
Reference  2.   These  approximations  are  listed  below. 
The  Standard  Normal  Cumulative  Distribution  Function 
Approximation : 

P(z)  =  1  -  s-t-(b,  +  t-(b2  +  t-(b,  +  t-(b4  +  t-b5)))) 
where   s  =  (1/ (2 -tt)  1/2)  -exp(-z2/2)   and  t  =  1/(1  +  b0-z).   And 
where 

b0  =  .2316419,   b,  =  .319381530,   b2  =  -.356563782, 

b4  =  -1.821255978,   b5  =  1.330274429, 
and   z  >  0.   For   z  <  0,  P(z)  =  1  -  P(|z|). 

The  Inverse  Standard  Normal  Cumulative  Distribution  Function 
Approximation : 

z(P)  =  t  -  (c0  +  t-(c,  +  t-c2))/(l  +  t-(d,  +  t-(d2  +  t-d3))) 
where   t  =  (ln(l/Q2)  ) 1/2  and  Q  =  1  -  P.   And  where 

c0  =  2.515517,   c,  =  .802853,   c2  =  .010328, 

d,  =  1.432788,   d2  =  .189269,   d3  =  .001308, 
and   .5  <  P  <  1.   For   0  <  P  <  .5,   z(P)  =  -z(l  -  P) . 
The  Inverse  Chi-Square  Cumulative  Distribution  Function 
Approximation : 

v(P|M)  =  M-[l  -  2/(9-M)  +  z  •  (2/(9 -M)  )1/2]3 
where   P(z)  =  P(v|M).   In  the  program,  the  inverse  standard 
normal  cumulative  distribution  function  approximation  is  used  to 
determine   z. 
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The  Noncentral  Chi-Square  Cumulative  Distribution  Function 
Approximation : 

P(w|M,E   sK2/a2)    =   P(z) 
where      z   =      [2-w/(l   +  b)]1/2      -    [2-a/(l   +  b)    -    1]1/2     with 
a   =  M  +   Z   sK2/a2,      b   =    (2   sK2/a2)/(M  +   2   sK2/a2)      and   the   sum   index 
K=l,    2,     •••    ,    M.       In   the  program,    the   standard   normal 
cumulative  distribution   function  approximation   is  used  to 
determine      P(z) . 
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Appendix  3.   The  Magnetic  Signal 

The  encounter  models  are  defined  by  the  following 
conditions:   A  submarine  magnetic  anomaly  field  is  a  magnetic 
dipole  field  that  is  superimposed  on  an  earth  magnetic  field  that 
is  constant  over  an  encounter  region.   A  magnetometer  magnetic 
signal  value  is  the  magnitude  of  the  projection  of  a  dipole 
magnetic  field  on  the  earth  magnetic  field  at  the  location  of  the 
magnetometer.   The  basis  for  determining  magnetic  signal  values 
in  the  program  is  an  expression  that  can  be  developed  as 
follows:    In  the  right  handed  rectangular  coordinate  system  that 
is  shown  in  Figure  1,  a  magnetic  dipole  is  at  the  origin,  the  xy- 
plane  is  the  horizontal  plane  at  a  representative  point  in  an 
encounter  region,  the  positive  y-axis  is  in  the  direction  of 
magnetic  north,  the  positive  x-axis  is  in  the  direction  of 
magnetic  east  and  the  positive  z-axis  is  positive  upward.   In 
this  rectangular  coordinate  system,  the  constant  earth  magnetic 
field  can  be  expressed  by   HE  =  HE-(j-cos  *  +  k-sin  $)   where   $ 
is  a  magnetic  dip  angle  and  HE   is  a  magnetic  field  magnitude 
that  characterizes  the  earth  field  in  an  encounter  region. 

In  a  spherical  coordinate  system  with  the  origin  at  the 
magnetic  dipole  and  the  polar  axis  in  the  direction  of  the  dipole 
moment,   Hd  =  (c-p/r3)  -(2-r-cos  9  +  ©-sin  6)   is  the  magnetic 
field  of  the  dipole  at  a  point  whose  spherical  coordinates  are 
(r,r,6).  In  this  expression  p   is  the  magnitude  of  the  dipole 
moment  and   c   is  a  constant  whose  value  is  determined  by  the 
choice  of  units.   The  magnetic  signal  for  a  magnetometer  that  is 
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Figure  1.  A  unit  vector  in  the  direction  of  the  earth  field  and  a 
unit  vector  in  the  direction  of  the  dipole  field  of  a 
dipole  at  the  origin  are  shown  at  a  point  that  is  a 
distance   r   from  the  dipole.   The  unit  vector   h   is 
in  the  direction  of  the  horizontal  component  of  the 
dipole  moment. 
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described  by  the  encounter  models  is   Hs  =  Ee-K(i/liE     when  the 
magnetometer  is  at  the  point  (r,r,8)  and   Hs  can  be  expressed  in 
terms  of  the  rectangular  coordinate  system  as  follows:   First, 
let  r0  be  a  unit  vector  in  the  direction  of  the  spherical 
coordinate  system  polar  axis.   Then  the  magnetic  dipole  field 

Hd  =  (c-p/r3)  -(3-r-cos  6  -  r0) 
since   r  =  (i-x  +  j-y  +  k-z)/r  where  r  =  (x2  +  y2  +  z2)1/2  , 
0  =  [  (r0  x  r)  x  r]/sin  9  =  [(r«r0)  -r  -  (r*r)  -r0]/sin  6,   and  the 
dot  product   r«r0  =  cos  0.   The  unit  vector  r0  can  be  expressed 
in  the  rectangular  coordinates  by  noting  that   r0  =  p/p  and  then 
expressing  p   in  rectangular  coordinates.   To  do  this,  let  ft 
be  the  depression  angle  of   p   from  the  xy-plane  (the  horizontal 
plane)  with  ft   positive  downward  and  let  a      be  the  direction  of 
p   relative  to  magnetic  north.   Then  the  magnetic  dipole  moment 
p  =  p-(h-cos  ft  -  k-sin  ft)   in  terms  of  the  unit  vector  k  and 
the  unit  vector  h  =  i-sin  a   +  j -cos  ft  which  has  the  direction 
of  the  horizontal  component  of  p.   With  these  results,  the  unit 
vector  r0  =  i  •  (cos  ft -sin  a)    +   j  •  (cos  ft -cos  a)  -  k-sin  ft  and 
Hs  =  (c  "P/1"3)  '  [  (3/r)  -cos  9  •  (y  -cos  $  -  z  -sin  $) 

-  (cos  §-cos  ft -cos  a  +  sin  §-sin  ft)] 
where   cos  6  =  (1/r) -(x-cos  ft -sin  a  +  y-cos  ft -cos  a  -  z-sin  ft) 
since   cos  6  =  r*r0.   As  can  be  seen  from  this  expression,  for  a 
constant  dipole  moment  magnitude  and  direction  and  a  constant 
earth  field  magnitude  and  direction,  the  magnetic  signal  is  only 
a  function  of  the  rectangular  coordinates  of  the  location  of  the 
magnetometer  relative  to  the  dipole.   In  the  encounter  models, 
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both  of  these  conditions  are  satisfied.   However,  by  allowing  p, 
fi,  a,  Hs  and   $   to  vary,  the  expression  for  Hs   is  applicable 
to  more  general  encounter  models. 
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Appendix  4.   The  Anderson  Formulation 

In  the  encounter  models,  the  magnetic  signal   Hs  at  a 
sample  point  in  a  straight  line  encounter  can  be  represented  in  a 
form  described  by  Anderson  in  Reference  4.   For  convenience,  the 
Anderson  formulation  is  used  in  the  program  to  determine  values 
for  Hs.   It  can  be  developed  as  follows:   In  Figure  2,  the 
primed  rectangular  coordinate  system  is  superimposed  on  the 
rectangular  coordinate  system  of  Figure  1  so  that  the  origin  is 
coincident  with  the  origin  of  that  system.   A  magnetometer  is  in 
a  straight  line  encounter  with  a  magnetic  dipole  that  is  located 
at  the  origin  of  the  combined  system.   The  combined  system  moves 
with  the  magnetic  dipole  with  the  x'-axis  oriented  so  that  it  is 
parallel  to  and  in  the  direction  <p     of  magnetometer's  relative 
to  the  magnetic  dipole  and  the  z'-axis  oriented  so  that  it  is 
directed  toward  and  passes  through  the  CPA  on  that  track.   Let 
1,  m,  and   n   be  the  direction  cosines  of  the  dipole  moment  p 
and   1,,  m,   and  n,   be  the  direction  cosines  of  the  earth 
magnetic  field  HE  in  the  primed  system.   Then  the  unit  vector 
r0  =  i'-l  +  j'-n  +  k1 -m.   And,  for  points  on  the  relative  track, 
x'  =  s',  y'  =  0   and   z'  =  R  where  R   is  the  slant  range  of  the 
dipole  at  CPA  and  where   s'  is  the  algebraic  distance  of  the 
magnetometer  from  CPA  on  the  relative  track.  (It  is  negative  for 
points  before  CPA  and  positive  for  points  after  CPA.)   This 
implies  that  magnetometer  position  vector  in  the  moving 
coordinate  system  is   r  =  i'-Js'/r)  +  k' • (R/r) . 
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Figure  2.  The  dipole  is  at  the  origin  in  the  unprimed  and  primed 
coordinate  systems.   In  the  primed  coordinate  system, 
the  CPA  is  at   (0,  0,  R)   and,  for  a  time   t   during 
the  encounter,  the  magnetometer  is  at   [s'(t),  0,  R] . 
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From  Appendix   3 , 

Hd  =    (c-p/r3)  •  (3   r-cos  e   -   r0)  . 
With      r      and      r0      expressed   in  terms   of   the  primed  unit  vectors 
and     cos  e  =  r«r0  =   l-(s'/r)    +  n*(R/r),    this  becomes 

Hd  =    (c-p/r3)  •  [  (3/r2)  • (1 -s"    +  n-R) • (i' -s'    +  k'  -R) 
-    (i' -1   +j ' -m  +   k' -n) ] . 
Then,    since        He/He  =   i'    1,   +   j  •  -m,  +  k'  -n,      and     Hs  =  KE'K(i/HEf 
Hs  =    (c-p/r3)  -[(2-1-1,   -  m-m,   -  n-n,)  -(s1)2 

+   3  •  (n  -1,   +   l-nJ-s'-R  +    (2  -n-n,   -   1-1,   -  m-m,)  -R2]  . 
The   quantities 

Aj,=    2-1-1,    -   m-m,    -   n-n, 

A,    =   3  •  (n-1,   +    1  -n,) 

Aq  =   2-n-n,    -    1-1,    -  m-m,, 

are  the  Anderson  coefficients.   With   r  =   [(s')2  +  R2]1/2   and 
6  =  s'/Rf   Hs   can  now  be  expressed  as  follows: 

HS(B)  =  (c-p/R3)  -Z  hyFAZ) 
where  the  FJR)      =      BJ/(1  +  B2)5/2   for  J  =  0,  1,  2   are  the 
Anderson  functions.   This  is  the  Anderson  formulation. 

To  relate  the  Anderson  formulation  for  Hs  to  the 
formulation  for  Hs   in  Appendix  3,  first  note  that 


1  =  r0-i' 


m  =  r0  •  j  ■ 


n  =  r0-k' 


and 
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1,  =  (He/He)  -i' 
m,  =  (He/He)- j« 

n,  =  (HE/HE)-k«. 

Then  express   r0,  (He/He)   and  the  unit  vectors   i',  j'  and  k' 
in  terms  of  the  unit  vectors   i,  j   and  X  and  take  the 
indicated  dot  products.   From  Appendix  3, 

r0  =  i  •  (cos  n-sin  a)  +  j  •  (cos  n-cos  a)  -  k-sin  n 
and 

He/He  =  j -cos  $  +  k-sin  $. 
To  express   i',  j'  and  k1  in  terms  of  the  unit  vectors   i,  j 
and  k,  note  that  the  unprimed  coordinate  system  can  be 
transformed  to  the  primed  coordinate  system  by  two  rotations  that 
are  defined  as  follows:   First,  rotate  a  coordinate  system  that 
is  coincident  with  the  unprimed  coordinate  system  about  the 
z-axis  through  the  angle  (0  -  n/2)    with  positive  angles  clockwise 
(left  hand  rule)  so  that  its  x-axis  is  parallel  to  and  in  the 
direction  of  the  relative  track.   Next,  rotate  this  system  about 
its  x-axis  through  an  angle  8      with  positive  angles 
counterclockwise  (right  hand  rule)  so  that  the  positive  z-axis 
passes  through  the  CPA.   The  angle  6      is  related  to  the 
vertical  separation   z0  of  the  magnetometer  and  the  dipole  and 
the  algebraic  encounter  lateral  range   L  that  is  positive  if  the 
dipole  is  to  the  left  of  the  relative  track.   With  these  sign 
definitions:   L  =  R-sin  6      and   z0  =  R-cos  <5 .   After  the 
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rotation,  the  auxiliary  coordinate  system  is  coincident  with  the 
primed  coordinate  system. 

These  transformations  can  be  described  in  terms  of  matrix 
equations  as  follows:   Let   (x",y",z" )   be  the  coordinates  of  a 
point  in  the  coordinate  system  that  is  coincident  with  the 
auxiliary  coordinate  system  after  the  first  rotation.   Then  the 
transformation  from  the  unprimed  coordinate  system  to  this  double 
primed  coordinate  system  is  described  by  the  matrix  equation 


X" 

y" 

= 

z" 

sin  0   cos  0    0 

■cos  0    sin  0    0 

0       0       1 


x 

y 

z 


And  the  transformation  from  the  double  primed  coordinate  system 
to  the  primed  coordinate  system  is  described  by  the  matrix 
equation 


X' 

y 

= 

z1 

10         0 

0     cos  6  sin  6 

0    -sin  S  cos  6 


X" 

y" 

z" 

Taking  the  product  of  the  rotation  matrices  in  the  indicated 
order  yields  the  matrix  equation 


X' 

y 

= 

z' 

sin  0  cos  0  0 

-cos  6 -cos  0    cos  6 -sin  0    sin  6 
sin  6  -cos  0   -sin  S -sin  0    cos  S 


x 

y 

z 


28 


which  defines  the  transformation  from  the  unprimed  to  the  primed 
coordinate  system.   The  unprimed  vector  components  of  the  unit 
vector   i'  can  be  found  by  transforming  the  coordinates  (1,0,0) 
in  the  primed  system  to  their  corresponding  coordinates  in  the 
unprimed  system  with  the  inverse  of  this  matrix  and  then 
repeating  this  process  for  (0,1,0)  and  (0,0,1)  in  order  to  find 
the  unprimed  unit  vector  components  of  j'   and  k1.   However, 
since  the  inverse  transformation  matrix  is  the  transpose  of  this 
matrix,  the  elements  of  the  row  that  corresponds  to  a  primed  unit 
vector  are  the  magnitudes  of  the  unprimed  vectors  that  are  its 
components.   Consequently: 

i'  =   i-sin  0  +  j -cos  <p 

ji  =  -i -cos  •  <5 -cos  <p  +   j -cos  6 -sin  <p  +   k-sin  6 
k1  =   i-sin  6-cos  <p   -  j -sin  <5-sin  0  +   k-cos  S. 
Then,  taking  the  dot  products  between  r,  r0  and  these  three 
unit  vectors  as  indicate  above  gives: 

1  =  cos  ft -cos  (0  -  a) 

m  =  cos  S  -cos  ft-sin  (0  -   a)  -   sin  <S-sin  ft 
n  =  -sin  6  -cos  n-sin  (<p   -   a)    -   cos  «S-sin  ft 
and 

1,  =  cos  $  -cos  0 

m,  =  cos  S  -cos  $-sin  0  -  sin  6  -sin  § 

n,  =  -sin  6  -cos  §-sin  0  -  cos  S -sin  $. 

These  are  the  relations  that  are  used  in  the  program  to  determine 
values  for  the  Anderson  coefficients. 
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Appendix  5.   The  Encounter  Equations  of  Motion 

In  the  double  primed  coordinate  system  that  is  defined  in 
Appendix  4,  the  equations  of  motion  of  a  magnetometer  relative  to 
a  submarine  (dipole)  are: 

x" (t)  =  s' (t) 

Y"(t)  =  -L 

2»(t)  =  Z0 

where   L  is  the  algebraic  encounter  lateral  range  that  is 
defined  in  Appendix  4 ,  z0   is  the  vertical  separation  between  the 
magnetometer  and  the  submarine  and   s'(t)   is  the  distance  of  the 
magnetometer  from  the  CPA  on  the  relative  track.   With  w  the 
speed  of  the  magnetometer  relative  to  the  submarine  and  t  a 
relative  time  parameter,   s'(t)  =  w-t  .   These  equations  can  be 
considered  to  be  the  ones  used  in  the  program  to  describe  the 
motion  of  a  magnetometer  relative  to  a  submarine.   There,   t   is 
determined  by  t  =  [J  -  (M-l)/2]  -<St  where  the  index  J  =  1,  2, 
•••,  M  and  St,    a  time  step,  is  the  time  between  samples.   Note, 
when   t  =  0,  the  magnetometer  is  at  the  CPA. 

In  the  coordinate  system  of  Figure  1  in  a  straight  line 
encounter  as  defined  in  the  encounter  models,  the  equations  of 
motion  of  a  magnetometer  relative  to  a  submarine  can  be  written 
as  follows: 

x(t)  =  s'  (t)  -sin  0   +   L-cos  0 

y(t)  =  s1 (t) -cos  0   -   L-sin  0 

z(t)  =  z0, 


30 


since  the  transformation  from  the  double  primed  coordinate  system 
to  the  primed  coordinate  system  is  determined  by  the  matrix 
equation 


x 

y 

z 


sin  0  -cos  0  0 
cos  0  sin  0  0 
0        0       1 


X" 

y" 

z" 

x"  =  s'(t)   and   y"  =  -L.   The  above  equations  and  the  expression 

for  Hs   in  Appendix  3  could  have  been  used  in  the  proqram  to 
evaluate  the  magnetic  signal.   In  particular,  with  these 
equations  of  motion  and  with  the  two  relations   L  =  R-sin  S      and 
z0  =  R-cos  6,    the  expression  for  Hs   in  Appendix  3  can  be  written 
in  terms  of   n,  a,  $,  0,  6      and  R   so  that  it  is  identical  in 
appearance  to  the  Anderson  formulation  for  Hs   in  terms  of  these 
quantities.   The  definition  of  6      in  Appendix  4  in  terms  of  a 
counterclockwise  rotation  results  in  a  definition  of  the 
algebraic  lateral  range  that  is  consistent  with  some  that  have 
been  used  elsewhere. 

In  the  program,  the  relation  w  =  v  -  u   is  the  basis  for 
determining  the  relative  speed  w.   In  this  relation,  v   is  the 
velocity  of  the  magnetometer,  u   is  the  velocity  of  the  submarine 
(dipole)  and  w   is  the  velocity  of  the  magnetometer  relative  to 
the  submarine.   This  relation  implies  the  following  equations: 
wx  =  vx  -  ux,   wy  =  vy  -  uy   and  wz  =  v2  -  uz  where  the  coordinates 
x,  y  and   z   refer  to  a  fixed  coordinate  system  with  the  same 
orientation  as  that  of  Figure  1.   In  the  encounters  of  the 
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models,  vx  =  v-sin  a,    vy  =  v-cos  a      and  vz  =  0   where  a    is  the 
magnetic  course  and  v   is   speed  of  the  magnetometer.   And,  in 
addition,   ux  =  u-sin  B,  Uy  =  u-cos  6   and   u2  =  0   where   6   is 
the  magnetic  course  and  u   is  the  speed  of  the  submarine.   The 
relative  magnetic  course  <p     and  the  relative  speed  of  the 
magnetometer  are  defined  by  wx  =  w-sin  <p     and  wy  =  w-cos  0.   In 
the  program,  <p      and  w   are  determined  with  a  rectangular  to 
polar  conversion  routine  where  w  =  (  wx2  +  wy2)1/2  and  where  <p 
is  determined  by   sin^wyw)   and   cos"1(w/w). 
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Appendix  6.   The  Submarine  Magnetic  Dipole 

In  the  encounter  models,  the  magnitude  and  direction  of  a 
submarine's  dipole  moment  p   are  determined  by  first  determining 
its  components  in  the  rectangular  coordinate  system  of  Figure  1. 
In  that  coordinate  system, 

Px  =  (Plp  +  Pu)  'sin  fi  +  (Ptp  +  Pn)  -cos  B> 

Py  =  (Plp  +  Pu)  -cos  6  -  (p^  +  pj  -sin  B 

Pz  =  "(Pvp  +  Pv.) 
where   B   is  the  submarine's  magnetic  course,  p^,  p^  and  pw 
are  the  permanent  and   pu,  p^   and  p^   are  the  induced 
longitudinal,  transverse  and  vertical  magnetic  dipole  moments  of 
the  submarine.   These  relations  are  based  on  the  following  sign 
convention: 

pL  is  positive  when  pL   is  directed  from  stern  to  bow. 

pT  is  positive  when  pT   is  directed  from  port  to  starboard. 

pv  is  positive  when  pv   is  directed  downward. 
The  permanent  dipole  moments  are  input  parameters  in  the  program 
and  the  induced  dipole  moments  are  determined  in  a  way  that  is 
similar  to  one  that  is  described  in  Reference  1.   In  the 
encounter  models,  a  submarine  is  a  ferromagnetic  prolate 
ellipsoid  with  the  major  axis  the  submarine's  longitudinal  axis 
and  the  equal  minor  axes  the  submarine's  transverse  and  vertical 
axes.   And,  the  induced  dipole  moments  are: 
Pu  =  kL-HEL 
P^  =  K-Ket 

Pvi   =   ^V  *  ^EV 
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where  HEL/  H^  and  H^  are  the  vector  components  of  the  earth 
magnetic  field  in  the  rectangular  coordinate  system  defined  by 
the  ellipsoid  axes  and  the  submarine  magnetic  moment  sign 
convention.   The  magnitude  of  these  vector  components  are: 

HEL  =      HE-cos   $-cos   B 

Het  =   -HE-cos   $-sin  B 

hev  =  HE-sin  $. 
The  earth  magnetic  field  dip  angle   *   and  the  earth  magnetic 
field  magnitude  HE  can  each  be  chosen  to  be  input  parameters  in 
the  program  or  the  can  be  computed  by  the  program  as  described  in 
Appendix  7 .  (Values  for  $   and  HE  can  be  found  from  magnetic 
data  charts,  for  example,  see  References  5  and  References  6.) 
By  using  the  above  relations, 

pK  =  HE-cos  #  •  (k,_  -  kT)  -sin  B-cos  B  +  (p^-sin  B  +  p-^-cos  B) 
py  =  HE-cos  §-(kL-cos2B  +  kT-sin2B)  +  (p^-cos  B  -  p^-sin  B) 
p2  =  -(Hg-kv-sin  $  +  pVP) 
In  the  encounter  models,  the  values  of  the  permeability 
coefficients   kL,  kT  and  k^  are  related  to  submarine 
displacement  by  the  following  relations: 

K   =  fL-w 

kT  =  fT-W 

K   =  fv-w 
where  W   is  the  submarine  displacement  in  tons  and   fL,  fT  and 
fv  are  permeability  factors  that  are  determined  by  a  submarine's 
magnetic  characteristics.   In  the  program,  the  units  of  HE  are 
gamma,  the  units  of  p  are  in  oersted-centimeter3,  the  units  of 
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k  are  oersted-centimeter3/gamma  and  the  units  of   f   are 
oersted-cent  imeter3/gamma-ton.   If  the  units  of   HE  were   gamma, 
but  the  units  of  p  were  gamma-foot3,  then  the  units  of  k 
would  be  foot3  and  the  units  of   f  would  be  foot3/ton.   To 
convert  p   in  gamma-foot3  to  oersted-centimeter3  or  k   in  foot3 
to  oersted-centimeter3/gamma  or   f   in  foot3/ton  to  oersted- 
centimeter3/gamma-ton,  divide  by  3.53.  (The  program  default 
values  for   fL,  fT  and   fv  are  values  from  Reference  1  in 
foot3/ton  that  have  been  divided  by  3.53  to  give  values  in 
oersted-cent imeter3/gamma-ton.  )   To  convert   p   in  weber-meter  to 
oersted-centimeter3,  multiply  by   [  1/  (4  -ir)  ]  -1010  and  to  convert 
p   in  ampere-meter2  to  oersted-centimeter3,  multiply  by   1-103. 
Note,  HE   in  oersted  equals   10'5  times  HE   in  gamma,  HE  in  gamma 
equals   HE   in  nanotesla  and  HE   in  tesla  equals   HE   in 
weber/meter2,  so  HE   in  gamma  equals   109  times  HE   in 
weber/meter2  and  HE   in  oersted  equals   104  times  HE   in 
weber/meter2. 


35 


Appendix  7.   The  Earth  Magnetic  Field  and  Dip  Angle 

An  auxiliary  magnetic  field  model  is  described  in  this 
appendix.   The  model  is  the  basis  for  a  default  choice  for  either 
the  value  of  the  earth  magnetic  field  magnitude  parameter  HE   or 
the  earth  magnetic  field  dip  angle  parameter   $.   Relative  to 
encounter  model  accuracy,  the  default  values  should  be  adequate 
in  most  cases. 

In  the  model,  the  earth  magnetic  field  is  generated  by  a 
magnetic  dipole  that  is  located  at  the  earth's  center  and  the 
earth  is  a  nonmagnetic  sphere  of  radius   rE.   With   pE   the 
magnitude  of  the  dipole  moment,  the  magnitude  of  the  earth  field 
at  a  point  is 

HE  =  (HE-HE)1/2  =  (c-pe/r3)  -(3-cos2e  +  1)1/2 
where  6   is  the  polar  angle  of  the  point  in  a  spherical 
coordinate  system  and   c   is  a  constant  whose  value  is  determined 
by  the  choice  of  units.   The  dipole  moment  is  coincident  with  and 
in  the  direction  of  the  polar  axis  which  is  directed  toward  the 
earth's  southern  hemisphere.   In  this  coordinate  system,  at  any 
point  on  the  surface  of  the  earth: 

HE  =  HE0-  (3  -cos29  +  1)1/2 
where  HE0  is  the  value  of  HE  at  the  magnetic  equator  which  is 
defined  by  the  points  on  the  earth  where  e  =  90°.   In  terms  of 
the  dip  angle,  at  any  point  on  the  earth's  surface: 

HE  =  2  -HE0-  (3  -cos2$  +  1)"1/2. 
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This  expression  can  be  obtained  by  noting  that   $   can  be 
defined  in  terms  of  the  or  the   r  and  ©   components   of  HE   in 
the  spherical  coordinates  as  follows: 

sin  $  =  -HEr/HE  =  -2-HE0-(cos  6)/HE 
and 

cos  $  =  H^H,  =  HE0-(sin  6)/HE. 
Based  on  these  relations, 

[(sin  $)/2]2  +  cos2$  =  (HE0)2 
and 

HE  =  HE0- [  (sin2)/4  +  cos2$]"1/2  =  2  -HE0  •  (3  -cos2$  +  1)"1/2 
The  dip  angle  is  determined  from  a  magnetic  latitude  for  the 
encounter  region. 

The  magnetic  latitude  and  longitude  of  the  representative 
point  of  an  encounter  region  can  be  defined  in  terms  of  its 
geographic  latitude  and  longitude  by  using  the  following 
transformations:   First,  convert  the  latitude  and  longitude  of 
the  point  to  rectangular  coordinates  in  a  right-handed  coordinate 
system  whose  origin  is  at  the  center  of  a  spherical  earth,  whose 
positive  z-axis  passes  through  its  north  geographic  pole  and 
whose  positive  y-axis  passes  through  the  point  on  its  equator 
with  latitude  0°  and  longitude  0°.   Next,  rotate  a  coordinate 
system  that  is  coincident  with  this  system  in  a  clockwise 
direction  (left  hand  rule)  about  its  z-axis  so  that  its  positive 
y-axis  passes  through  the  point  with  latitude  0°  and  with 
longitude  equal  to  that  of  the  north  magnetic  pole.   Then,  rotate 
the  system  in  a  clockwise  direction  about  its  x-axis  so  that  its 
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positive  z-axis  passes  through  the  north  magnetic  pole.   Next, 
transform  the  rectangular  coordinates  of  the  representative  point 
in  this  system  to  its  coordinates  in  the  spherical  coordinate 
system  that  is  associated  with  it.   Then,  with  n  the  polar 
angle  of  the  representative  point  in  this  system,  the  magnetic 
latitude  of  the  representative  point  is  1^   =  90°  -  n.   Since 
tan  §  =  -2 -cot  6   from  above  and  9  =  180°  -  n,  the  dip  angle   * 
is  given  by  the  following  relation: 

$  =  tan'1  (2  -tan  1^)  . 
In  the  program,  the  transformations  described  above  are 
accomplished  in  part  by  a  rectangular  to  polar  conversion 
routine.   In  particular,  by  rotating  the  final  rectangular 
coordinate  system  about  its  z-axis  so  that  the  x-coordinate  is 
zero,  the  polar  angle  of  the  polar  coordinates  of  the 
representative  point  in  the  resulting  yz-plane  determines   L^. 

Some  values  of   $   and   HE   that  are  listed  in  Table  1  were 
generated  using  the  program.   The  latitude  and  longitude  of  the 
magnetic  pole  that  are  in  the  program  are   76°  N  and  100°  W. 
Values  of   $   from  Reference  6  and  of  HE   from  Reference  7  are 
also  listed  in  Table  1.   A  comparison  of  the  program  values  with 
these  values  gives  an  indication  of  the  errors  inherent  in  the 
procedure. 
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Encounter 

Encounter 

Program 

Chart  #30 

Program 

Chart  #39 

Latitude 

Longitude 

Dip  $ 

Dip  § 

Field  HE 

Field  HE 

60°N 

180°W 

74° 

70° 

.63 

.53 

60°N 

150°W 

78° 

73° 

.66 

.55 

60°N 

90°W 

82° 

84° 

.68 

.61 

60°N 

30°W 

75° 

74" 

.64 

.52 

30°N 

150°W 

58° 

50° 

.51 

.41 

30°N 

90°W 

62° 

60° 

.55 

.51 

30°N 

60°W 

59° 

59° 

.53 

.48 

30°N 

30°E 

37° 

43° 

.41 

.42 

30°N 

60°E 

31° 

45° 

.39 

.46 

30°N 

90°E 

30° 

44° 

.39 

.44 

308N 

150°E 

42° 

40° 

.43 

.40 

0° 

90°W 

26° 

18° 

.38 

.32 

0° 

30°W 

9° 

-5° 

.35 

.28 

0° 

0° 

-48° 

-25° 

.35 

.31 

0° 

60°E 

-25° 

-20° 

.38 

.37 

30°S 

180°W 

-45° 

-55° 

.44 

.50 

30°S 

90°W 

-30° 

-34° 

.39 

.31 

30°S 

90°E 

-62° 

-66° 

.55 

.52 

60°S 

30°W 

-69° 

-56° 

.60 

.34 

Table  1.   Some  corresponding  values  for  the  encounter  dip  angle 
$   and  the  earth  magnetic  field  intensity  HE   in 
oersted.   Chart  #33  refers  to  Reference  5  and  Chart  #39 
refers  to  Reference  6. 

A  small  difference  between  a  chart  value  and  a  computed  value  for 
the  dip  angle   $   at  some  location  does  not  imply  a  small 
difference  between  the  chart  value  and  the  computed  value  for  the 
earth  magnetic  field  intensity  at  that  location.   For  example, 
note  the  values  for   60°N   latitude  and   30°W   longitude. 
(Positive  values  of   $   indicate  the  inclination  or  dip  below  the 
horizontal  of  the  north  seeking  end  of  a  dip  needle.   Negative 
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values  indicate  the  inclination  below  the  horizontal  of  the  south 
seeking  end.) 

An  extension  of  the  above  procedure  can  be  made  for  finding 
the  magnetic  variation  at  a  location.   However,  the  values 
generated  by  using  the  procedure  are  generally  unsatisfactory. 
Magnetic  variation  values  are  charted  in  Reference  7. 
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Appendix  8.   An  Alternative  Encounter  Model 

The  alternative  model  that  is  referred  to  in  Section  III  is 
described  in  more  detail  in  this  appendix.   In  the  model,  the 
detection  range  is  the  slant  range   R  at  the  CPA  for  an 
encounter  with  a  specified  detection  probability  (usually  .5)  is 
defined  by: 

R  =  [c-p/Hs],/3 
where   c   and   p   are  defined  in  Appendix  3   and  Hs  represents 
a  minimum  detectable  average  magnetic  signal  that  is  defined  by: 

Hs  =  (ORF) -NM 
where   ORF  is  a  signal-to-noise  ratio  called  the  operator 
recognition  factor  and  NM   represents  the  magnetic  noise. 
Combining  these  two  relations  gives: 

R  =  <c-p/[(ORF)  -NM]}1/3. 
The  value  for   ORF   depends  both  on  the  specified  detection 
probability  and  on  a  specified  or  implied  false  alarm 
probability. 

The  Anderson  formulation  is  consistent  with  these  relations 
in  an  approximate  sense  if  the  average  magnetic  signal   H   is 
defined  as  a  root  mean  square  value  such  that  H  =  (c-p/R3)  -k 
where   k   is  an  encounter  parameter  defined  by: 

k  =  (2  [S  A.-F^BJ]2}172 
with  the  first  sum  index   K  =  1,  2,  •••  ,  M   and  the  second  sum 
index  J  =  0,  1,  2.   For  a  particular  encounter  geometry,  k   is 
constant  and  this  suggests  that  the  two  encounter  models  could  be 
used  to  determine  an  average  value  for   k   for  an  encounter 
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region  based  on  average  submarine  magnetic  characteristics. 
Values  for  both  k  and  H  are  generated  by  the  program  and  such 
values  can  give  an  indication  of  the  magnitude  of  the  differences 
in  detection  range  estimates  that  are  based  on  this  model  and 
those  that  are  based  on  either  of  the  other  two  encounter  models. 
A  more  detailed  comparison  of  these  encounter  models  is  described 
in  Reference  8. 
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Appendix  9.   An  Example  of  the  Program  Output 

The  program  is  designed  to  generate  values  for  the  following 
quantities:  encounter  parameters,  lateral  range  functions  for  the 
crosscorrelation  encounter  model  and  for  the  square  law  encounter 
model,  average  magnetic  signal,  slant  range  at  CPA,  encounter 
parameter  values  and  magnetic  signal  values.   The  values  can  be 
saved  as  a  program  data  file  and/or  they  can  be  printed.   Tables 
2  through  7  are  examples  of  the  program's  printed  output. 

Figures  3  and  5  are  plots  of  the  lateral  range  function 
values  that  are  listed  in  Tables  3  and  6  for  the  square  law 
encounter  model.   For  either  the  square  law  or  the 
crosscorrelation  encounter  model,  as  the  magnitude  of  the  lateral 
range  increases,  lateral  range  function  values  (detection 
probabilities)  do  not  approach  zero  but  approach  the  value  of  the 
encounter  false  alarm  probability.   A  necessary  condition  for  the 
formal  definition  of  sweep  width  given  in  Reference  9  to  be 
meaningful  is  that  the  lateral  range  function  approach  zero  as 
the  lateral  range  increases  without  limit.   Since  this  is  not  the 
case  for  either  encounter  model,  the  formal  definition  of  sweep 
width  must  be  modified  with  their  use.   Reference  10  provides  an 
example  of  this. 

Figures  4  and  6  are  plots  of  the  magnetic  signal  values  that 
are  listed  in  Tables  4  and  7. 
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program  file  name 

program  data  file  name 

magnetic  data  file  name 

processing  data  file  name 

kinematic  data  file  name 

combined  magnetic,  processing  &  kinematic  data  file  name  data.mpk 

encounter  latitude  (decimal  degrees) 

encounter  longitude  (decimal  degrees) 

encounter  variation  (decimal  degrees) 

encounter  dip  angle  (decimal  degrees) 

encounter  magnetic  field  intensity  (gamma) 

permanent  longitudinal  moment  (oersted-cm3) 

permanent  transverse  moment  (oersted-cm3) 

permanent  vertical  moment  (oersted-cm3) 

target  displacement  (tons) 

target  longitudinal  permeability  coefficient 

target  transverse  permeability  coefficient 

target  vertical  permeability  coefficient 

target  longitudinal  permeability  factor 

target  transverse  permeability  factor 

target  vertical  permeability  factor 

sampling  period  (seconds) 

integration  time  (seconds) 

adjusted  integration  time  (seconds) 

number  of  samples  per  encounter 

magnetometer  course  (decimal  degrees) 

magnetometer  speed  (knots) 

magnetometer  altitude  (meters) 

target  course  (decimal  degrees) 

target  speed  (knots) 

target  depth  (meters) 

magnetometer  relative  course  (decimal  degrees) 

magnetometer  relative  speed  (knots) 

magnetometer-target  vertical  separation  (meters) 

target  induced  longitudinal  dipole  moment  (oersted-cm3) 

target  induced  transverse  dipole  moment  (oersted-cm3) 

target  induced  vertical  dipole  moment  (oersted-cm3) 

magnetic  dipole  moment  (oersted-cm3) 

dipole  moment  azimuth  (decimal  degrees) 

dipole  moment  depression  angle  (decimal  degrees) 

distance  between  samples  on  the  relative  track  (meters) 

false  alarm  rate  (false  alarms  per  hour) 

false  alarm  probability 

magnetic  noise  (gamma) 

maximum  lateral  range  (meters) 

lateral  range  step  (meters) 

number  of  lateral  range  function  values 


ru\u .  at\i> 

data. mad 

data. mag 

data. pro 

data. kin 

data.mpk 

45 

-45 

-25 

68.41608 

59035.05 

0 

0 

0 

1800 

13140 

2880 

2880 

7.3 

1.6 

1.6 

.5 

20 

20.5 

41 

290 

180 

100 

20 

10 

100 

286.8202 

180.2776 

200 

2.017795E+08 

-4.422564E+07 

1.58099E+08 

2.601272E+08 

7.637512 

37.42885 

46.3714 

3 

1.666667E-02 

.35 

500 

20 

51 

Table  2.   An  example  of  an  encounter  parameter  values  printout, 
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L 

p(cc) 

P(sl) 

H 

R 

V- 

meters 

gamma 

meters 

-500 

.2155959 

.0255824 

.184247 

538.5165 

2  .818488 

-480 

.2721803 

. 029692 

.21034 

520 

2.879916 

-460 

. 3464673 

3.582633E- 

02 

.2410144 

501.5974 

2.943744 

-440 

.4415205 

4 .529992E- 

02 

.2785504 

483.3218 

3.009938 

-420 

.5573696 

6. 047795E- 

02 

.3231911 

465. 1881 

3  .078418 

-400 

.6872897 

.0856887 

.3764825 

447.2136 

3  .  149025 

-380 

. 8146549 

. 1287804 

.4404623 

429.4182 

3.221497 

-360 

.9159235 

.2030898 

.5177304 

411.8252 

3.29544 

-340 

.9745625 

.3272475 

.6108596 

394 .4617 

3.370274 

-320 

.9958122 

.5141411 

.7233319 

377.3593 

3.445185 

-300 

.9997211 

.7388175 

.8597726 

360.5551 

3.519032 

-280 

.9999951 

.9177597 

1.025158 

344 . 093 

3.590261 

-260 

1 

.9896971 

1.224529 

328.0244 

3.  656775 

-240 

1 

.999714 

1.463644 

312.41 

3.71578 

-220 

1 

.9999992 

1.74981 

297.3214 

3.763606 

-200 

1 

1 

2.085394 

282.8427 

3.795517 

-180 

1 

1 

2.470496 

269.0725 

3  .805514 

-160 

1 

1 

2.901352 

256.125 

3.786234 

-140 

1 

1 

3.359555 

244.1311 

3.72898 

-120 

1 

1 

3.819509 

233.2381 

3.624058 

-100 

1 

1 

4.239317 

223.6068 

3 .461596 

-80 

i 

1 

4 .532758 

215.4066 

3 . 232958 

-60 

1 

1 

4 .618943 

208.8061 

2.932828 

-40 

1 

1 

4  .44342 

203  .9608 

2.5617 

-20 

1 

1 

3.97739 

200.9975 

2. 128285 

0 

1 

1 

3 . 332139 

200 

1.651697 

20 

1 

.9999993 

2 .540665 

200.9975 

1. 165581 

40 

.  9S99931 

.908713 

1.76834 

203.9608 

.7394171 

60 

.9933111 

. 467368 

1. 110188 

208.8061 

.5637183 

80 

.  9998733 

.7874625 

.9039734 

215.4066 

.  7782836 

100 

.9999999 

.9856169 

1.086971 

223.6068 

1. 134662 

120 

1 

.9992399 

1.254863 

233.2381 

1. 490748 

140 

1 

.9998478 

1.308583 

244  .  1311 

1.811734 

160 

1 

.9998431 

1.29482 

256. 125 

2 . 089946 

180 

1 

.9994995 

1.230946 

269.0725 

2.326072 

200 

1 

.9974561 

1. 131703 

282.8427 

2.523877 

220 

1 

.9872744 

1.017684 

297.3214 

2.688154 

240 

.999999 

.9510013 

.9018685 

312.41 

2.823785 

260 

.9999747 

.8646194 

.7917129 

328.0244 

2.935307 

280 

.9996544 

.7242736 

.6909791 

344.093 

3. 026728 

300 

.9973592 

.5582789 

.6011183 

360.5551 

3.10149 

320 

.9876337 

.4037092 

.5222269 

377.3593 

3. 162495 

340 

.9608865 

.2816961 

.4536729 

394.4617 

3.212155 

360 

.9086407 

.1948741 

.3944811 

411.8252 

3.252463 

380 

.8301703 

.1365131 

.3443362 

429.4182 

3.285058 

400 

.7333605 

9.820431E- 

-02 

.3014853 

447.2136 

3  .  311281 

420 

.6299387 

.  0731037 

.2645045 

465.1881 

3  .332228 

440 

.  53C2704 

5.646636E- 

-02 

.2325951 

483.3218 

3.348799 

460 

.4408869 

4 . 522933E- 

-02 

.205043 

501.5974 

3.361725 

480 

.  3645446 

3 .747106E' 

-02 

. 1812233 

520 

3 . 371603 

500 

.3013791 

3.  199235E- 

-02 

.1605954 

538.5165 

3  .  378925 

Table  3.   An  example  of  a  lateral  range  function  values  printout 
The  heading  for  the  crosscorrelation  values  is  p(cc) 
and  the  heading  for  the  square  law  values  is  p(sl). 
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Figure  3.  A  plot  of  the  square  law  lateral  range  function  values 
that  are  listed  in  Table  3.   The  horizontal  axis  is 
encounter  horizontal  range  at  CPA  in  meters.   The 
vertical  axis  is  encounter  detection  probability. 
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data. mad  magnetic  signal  values 
for  a  lateral  range  of   40   meters 

relative  CPA  distance  magnetic  signal 

s  in  meters  Hs  in  gamma 

-927.4282  -1 . 619177E-02 

-881.0568  -1.834315E-02 

-834.6854  -2 . 08614 4E-02 

-788.314  -2.381637E-02 

-741.9426  -2.728804E-02 

-695.5712  -.0313641 

-649.1998  -3.613118E-02 

-602.8283  -4.165412E-02 

-556.4569  -4 . 792992E-02 

-510.0855  -5.478963E-02 

-463.7141  -6.169327E-02 

-417.3427  -6.730653E-02 

-370.9713  -.0686387 

-324.5999  -5 . 933568E-02 

-278.2285  -2 . 6523 14E-02 
-231.8571  .0540608 

-185.4856  .2215785 

"139.1142  .5163918 

-92.74282  .9104893 

-46.37141  1.184586 

0  1.003207 

46.37141  .3837067 

92.74282  -.2252113 

139.1142  -.5323231 

185.4856  -.5837545 

231.8571  -.5168925 

278.2285  -.4202569 

324.5999  -.3304064 

370.9713  -.2572838 

417.3427  -.2007385 

463.7141  -.1577988 

510.0855  -.125285 

556.4569  -.1005532 

602.8283  -8.158605E-02 

649.1998  -6.689598E-02 

695.5712  -5.539922E-02 

741.9426  -4. 630719E-02 

788.314  -3.904362E-02 

834.6854  -3  .  3 18428E-02 

881.0568  -.0284142 

927.4282  -2 . 449729E-02 


Table  4.   An  example  of  a  magnetic  signal  values  printout 
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1  4 

1.3  - 

1.2  - 

1.1  - 

1 

0.9  - 

0.8  - 

0.7  - 

0.6  - 

0.5  - 
0  4- 

0.3  - 

0.2  - 

0.1  - 


-0  1 
-0.2  - 
-0  3  - 
—0.4 

-0  5  - 
-0  6  - 
-0.7 


°  ^^^W-EK^HB-^^e 


-i — i — i — i — i — r 


-927.43 


-I — i— | — i — 
-463.71 


i  i — i — i — i — i — i — i  i 


463.71 


927.43 


Figure  4.  A  plot  of  the  magnetic  signal  values  listed  in  Table  4 
The  horizontal  axis  is  relative  CPA  distance  in  meters 
The  vertical  axis  is  magnetic  signal  in  gamma. 
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program  file  name  MAD.BAS 

program  data  file  name  datal.mad 

magnetic  data  file  name  data. mag 

processing  data  file  name  data. pro 

kinematic  data  file  name  datal.kin 
combined  magnetic,  processing  &  kinematic  data  file  name  datal.mpk 

encounter  latitude  (decimal  degrees)  45 

encounter  longitude  (decimal  degrees)  -45 

encounter  variation  (decimal  degrees)  -25 

encounter  dip  angle  (decimal  degrees)  68.41608 

encounter  magnetic  field  intensity  (gamma)  59035.05 

permanent  longitudinal  moment  (oersted-cm3)  0 

permanent  transverse  moment  (oersted-cm3)  0 

permanent  vertical  moment  (oersted-cm3)  0 

target  displacement  (tons)  1800 

target  longitudinal  permeability  coefficient  13140 

target  transverse  permeability  coefficient  2880 

target  vertical  permeability  coefficient  2880 

target  longitudinal  permeability  factor  7.3 

target  transverse  permeability  factor  1.6 

target  vertical  permeability  factor  1.6 

sampling  period  (seconds)  .5 

integration  time  (seconds)  20 

adjusted  integration  time  (seconds)  20.5 

number  of  samples  per  encounter  41 

magnetometer  course  (decimal  degrees)  290 

magnetometer  speed  (knots)  180 

magnetometer  altitude  (meters)  100 

target  course  (decimal  degrees)  273 

target  speed  (knots)  10 

target  depth  (meters)  100 

magnetometer  relative  course  (decimal  degrees)  290.9828 

magnetometer  relative  speed  (knots)  170.462 

magnetometer-target  vertical  separation  (meters)  200 
target  induced  longitudinal  dipole  moment  (oersted-cm3)    1.339682E+08 

target  induced  transverse  dipole  moment  (oersted-cm3)  5. 522353E+07 

target  induced  vertical  dipole  moment  (oersted-cm3)  1.58099E+08 

magnetic  dipole  moment  (oersted-cm3)  2.144584E+08 

dipole  moment  azimuth  (decimal  degrees)  295.4022 

dipole  moment  depression  angle  (decimal  degrees)  47.49356 
distance  between  samples  on  the  relative  track  (meters)    43.84662 

false  alarm  rate  (false  alarms  per  hour)  3 

false  alarm  probability  1.666667E-02 

magnetic  noise  (gamma)  .35 

maximum  lateral  range  (meters)  500 

lateral  range  step  (meters)  20 

number  of  lateral  range  function  values  51 


Table  5.   An  example  of  an  encounter  parameter  values  printout 
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datal.mad  lateral  range  function  values 


L 

p(cc) 

p(sl) 

H 

R 

k 

meters 

ganuna 

meters 

-500 

9.718663E-02 

1.815607E- 

-02 

.1162113 

538 

.5165 

2 

.  117317 

-480 

.1097076 

1.887764E- 

-02 

.1307452 

520 

2 

.066198 

-460 

. 1252767 

1.979502E- 

-02 

.1476683 

501 

.5974 

2 

.017377 

-440 

.1449347 

2.098504E- 

-02 

.1677343 

483 

.3218 

1 

.971743 

-420 

.1701671 

2.256657E- 

-02 

.1915806 

465 

.1881 

1 

.930378 

-400 

.2031047 

2.473001E- 

-02 

.2196898 

447 

.2136 

1 

.89458 

-380 

.2467793 

2.779396E- 

-02 

.2543952 

429 

.4182 

1 

.865873 

-360 

.3053632 

3.231674E- 

-02 

.2951223 

411 

.8252 

1 

.845999 

-340 

.3841007 

3.932704E- 

-02 

.3461031 

394 

.4617 

1 

.836884 

-320 

.4881344 

5.082678E- 

-02 

.4085031 

377 

.3593 

1 

.840555 

-300 

.6185669 

7.091362E- 

-02 

.4830141 

360 

.5551 

1 

.859016 

-280 

.7643011 

.1082717 

.5793262 

344 

.093 

1 

.89408 

-260 

.8946846 

.1810074 

.6993423 

328 

.0244 

1 

.947158 

-240 

.9731281 

.3214172 

.8450639 

312, 

.41 

2 

.019043 

-220 

.9973629 

.5584107 

1.022772 

297, 

.3214 

2 

.109683 

-200 

.9999472 

.8329318 

1.262045 

282 

.8427 

2 

.217957 

-180 

.9999999 

.9794214 

1.555326 

269 

.0725 

2 

.341425 

-160 

1 

.9996887 

1.906387 

256, 

.125 

2, 

.476081 

-140 

1 

.9999999 

2.315586 

244, 

.1311 

2, 

.616057 

-120 

1 

1 

2.774995 

233, 

.2381 

2, 

.753399 

-100 

1 

1 

3.289594 

223, 

.6068 

2 

.878014 

-80 

1 

1 

3.848771 

215, 

,4066 

2, 

.978009 

-60 

1 

1 

4.356607 

208, 

.8061 

3, 

,040665 

-40 

1 

1 

4.738054 

203. 

.9608 

3, 

.054164 

-20 

1 

1 

4.921432 

200. 

,9975 

3, 

.009893 

0 

1 

1 

4.862248 

200 

2. 

,904766 

20 

1 

1 

4.561756 

200, 

.9975 

2. 

,742697 

40 

1 

1 

4.06876 

203. 

.9608 

2, 

.534609 

60 

1 

1 

3.462592 

208. 

,8061 

2. 

,296925 

80 

1 

1 

2.876776 

215. 

,4066 

2, 

,049263 

100 

1 

.9999861 

2.341757 

223. 

,6068 

1. 

,812213 

120 

1 

.9908414 

1.862014 

233. 

,2381 

1, 

,60562 

140 

.9999632 

.8492265 

1.463262 

244. 

,1311 

1, 

,446881 

160 

.9973499 

.55795 

1.150695 

256. 

,125 

1, 

,348197 

180 

.9773154 

.3394268 

.9187382 

269. 

.0725 

1, 

.31296 

200 

.9311967 

.2232209 

.7327372 

282. 

,8427 

1, 

,334314 

220 

.8711751 

.1618173 

.5979505 

297. 

,3214 

1. 

.39856 

240 

.8071772 

.1253278 

.4948332 

312. 

,41 

1. 

,49089 

260 

.7414287 

.1006817 

.4163041 

328. 

,0244 

1. 

,599193 

280 

.6739562 

8.250537E- 

■02 

.3566572 

344. 

,093 

1. 

,714996 

300 

.6054811 

.0684997 

.3130742 

360. 

5551 

1. 

,832886 

320 

.5376865 

5.753178E- 

■02 

.2806611 

377. 

3593 

1. 

949611 

340 

.4726234 

4.891159E- 

•02 

.2530364 

394. 

4617 

2. 

,063321 

360 

.4121077 

4.213673E- 

-02 

.2280127 

411. 

8252 

2. 

173046 

380 

.3573765 

3.681085E- 

■02 

.2054994 

429. 

4182 

2. 

27835 

400 

.3090144 

3.261634E- 

•02 

. 1853368 

447. 

2136 

2. 

379115 

420 

.2670518 

.0293013 

.1673323 

465. 

1881 

2. 

475412 

440 

.2311289 

2.666843E- 

•02 

.151282 

483. 

3218 

2. 

567413 

460 

.2006615 

2.456539E- 

■02 

.1369845 

501. 

5974 

2. 

655341 

480 

.1749714 

.0228749 

.12437 

520 

2. 

739442 

500 

.1533755 

2. 150718E- 

■02 

.1132715 

538. 

5165 

2. 

819965 

Table  6.   An  example  of  a  lateral  range  function  values  printout 
The  heading  for  the  crosscorrelation  values  is  p(cc) 
and  the  heading  for  the  square  law  values  is  p(sl). 
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Figure  5.   A  plot  of  the  square  law  lateral  range  function  values 
that  are  listed  in  Table  6.   The  horizontal  axis  is 
encounter  horizontal  range  at  CPA  in  meters.   The 
vertical  axis  is  encounter  detection  probability. 
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datal.mad  magnetic  signal  values 
for  a  lateral  range  of   40   meters 

relative  CPA  distance 
s  in  meters 

-876.9324 

-833.0858 

-789.2391 

-745.3925 

-701.5459 

-657.6993 

-613.8527 

-570.006 

-526.1594 

-482.3128 

-438.4662 

-394.6196 

-350.7729 

-306.9263 

-263.0797 

-219.2331 

-175.3865 

-131.5399 

-87.69324 

-43.84662 

0 

43.84662 

87.69324 

131.5399 

175.3865 

219.2331 

263.0797 

306.9263 

350.7729 

394.6196 

438.4662 

482.3128 

526.1594 

570.006 

613.8527 

657.6993 

701.5459 

745.3925 

789.2391 

833.0858 

876.9324 


magnetic  si 

gnal 

Hs  in  gamma 

6.869243E- 

03 

8.988515E- 

03 

1.181241E- 

02 

1.561336E- 

02 

2.078687E- 

02 

2.791575E- 

•02 

3.787262E- 

•02 

5.198577E- 

■02 

7.231193E- 

■02 

.1020894 

.1464967 

.2139163 

.317984 

.480659 

.735856 

1.130141 

1.707647 

2.448238 

3.132554 

3.263226 

2.439828 

1.023829 

-.1487533 

-.7045798 

-.8055338 

-.7078842 

-.5631492 

-.4307948 

-.3256765 

-.2465582 

-.1881302 

-.1451051 

-.1132614 

-8.948089E- 

-02 

-7.152954E- 

-02 

-5.782279E- 

-02 

-4.723642E- 

-02 

-3.896839E' 

-02 

-3.244173E' 

-02 

-2.723733E- 

-02 

-2.304761E 

-02 

Table  7.   An  example  of  a  magnetic  signal  values  printout 
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Figure  6.   A  plot  of  the  magnetic  signal  values  listed  in  Table 
7.   The  horizontal  axis  is  relative  CPA  distance  in 
meters.   The  vertical  axis  is  magnetic  signal  in 
gamma . 
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Appendix  10.  A  Program  Listing 


10  CLS  :  CLEAR 

20  PRINT  "Magnetic  Anomaly  Detection  (MAD)  Lateral  Range  Function  Program" 

30  DIM  X0(70),  R0(70),  H(70),  HS(70,  150),  K(70),  PDCC(70),  PDSL(70) 

40  PI  =  3.141592654#:  CON  =  2  *  PI  /  360:  KON  =  1852  /  3600:  N$  =  "MAD.BAS" 

50  Q0  =  .2316419:  Q1  =  .31938153#:  Q2  =  - .356563782#:  Q3  =  1.781477937* 

60  04  =  -1.821255978#:  Q5  =  1.330274429# 

70  11  =  2.515517:  12  =  .802853:  13  =  .010328:  14  =  1.432788:  15  =  .189269:  16  =  .001308 

80  PRINT 

90  INPUT  "generate  data  or  print  a  program  data  file  (g/p)";  A$ 

100  IF  AS  =  "G"  OR  A$  =  "g"  THEN  GOTO  120 

110  IF  A*  ■  "P"  OR  AS  =  "p"  THEN  GOTO  3070  ELSE  GOTO  90 

120  PRINT  :  AS  =  "a" 

130  INPUT  "magnetic,  processing  &  kinematic  data  entry  by  combined  file  (y/n)";  AS 

140  IF  AS  =  "Y"  OR  A$  =  "y"  THEN  GOTO  2120 

150  IF  AS  =  "N"  OR  AS  =  "n"  THEN  GOTO  160  ELSE  GOTO  130 

160  INPUT  "magnetic  data  entry  by  file/keyboard  (f/k)";  AS 

170  IF  AS  =  "F"  OR  AS  =  "f"  THEN  GOTO  780 

180  IF  AS  =  "K"  OR  AS  =  "k"  THEN  GOTO  190  ELSE  GOTO  160 

190  INPUT  "latitude  in  decimal  degrees  (N  +)";  LAT:  LATR  =  LAT  *  CON 

200  INPUT  "longitude  in  decimal  degrees  (E  +)";  LNG:  LNGR  =  -LNG  *  CON 

210  INPUT  "variation  in  decimal  degrees  (E  +)";  DEC:  DECR  =  DEC  *  CON 

220  AS  =  "a" 

230  INPUT  "input  dip  angle  (y/n)";  AS 

240  IF  AS  =  "Y"  OR  AS  =  "y"  THEN  GOTO  260 

250  IF  AS  =  "N"  OR  AS  =  "n"  THEN  GOTO  290  ELSE  GOTO  230 

260  INPUT  "dip  angle  in  decimal  degrees  (north  magnetic  hemisphere  +)";  DIP 

270  DIPR  =  DIP  *  CON 

280  GOTO  390 

290  LATPR  =  76  *  CON:  LNGPR  =  100  *  CON 

300  D  =  SIN(LNGR  -  LNGPR)  *  COS(LATR):  E  =  COS(LNGR  -  LNGPR)  *  COS(LATR):  F  =  SIN(LATR) 

310  X  =  E:  Y  =  F 

320  GOSUB  4170 

330  T  =  T  -  (90  *  CON  -  LATPR):  E  =  R  *  SIN(T):  F  =  R  *  COS(T) 

340  X  =  E:  Y  =  D 

350  GOSUB  4170 

360  X  =  F:  Y  =  R 

370  GOSUB  4170 

380  DIPR  =  ATN(2  *  TAN(T)):  DIP  =  DIPR  /  CON 

390  AS  =  "a" 

400  INPUT  "input  encounter  magnetic  field  intensity  (y/n)";  AS 

410  IF  AS  =  "Y"  OR  AS  =  "y"  THEN  GOTO  430 

420  IF  AS  =  "N"  OR  AS  =  "n"  THEN  GOTO  450  ELSE  GOTO  400 

430  INPUT  "encounter  magnetic  field  intensity  in  gamma";  HE 

440  GOTO  460 

450  HE  =  70000!  /  SQR(3  *  COS(DIPR)  *  COS(DIPR)  +  1) 

460  PLM  =  0:  PTM  =  0:  PVM  =  0 

470  AS  =  "a" 

480  INPUT  "input  target  permanent  dipole  moments  (y/n)";  AS 

490  IF  AS  =  "Y"  OR  AS  =  "y"  THEN  GOTO  510 

500  IF  AS  =  "N"  OR  AS  =  "n"  THEN  GOTO  540  ELSE  GOTO  480 

510  INPUT  "permanent  longitudinal  moment  in  oersted-cm3  (stern-to-bow  +)";  PLM  520  INPUT  "permanent 

transverse  moment  in  oersted-cm3  (port-to-starboard  +)";  PTM 

530  INPUT  "permanent  vertical  moment  in  oersted-cm3  (downward  +)";  PVM 

540  INPUT  "target  displacement  in  tons";  WT 

550  AS  =  "a" 

560  INPUT  "input  target  permeability  coefficients  or  factors  (c/f)";  AS 

570  IF  AS  =  "C"  OR  AS  =  "c"  THEN  GOTO  590 

580  IF  AS  =  "F"  OR  AS  =  "f"  THEN  GOTO  630  ELSE  GOTO  560 

590  INPUT  "longitudinal  permeability  coefficient  in  cgs  units";  KL 

600  INPUT  "transverse  permeability  coefficient  in  cgs  units";  KT 

610  INPUT  "vertical  permeability  coefficient  in  cgs  units";  KV 

620  FL  =  KL  /  WT:  FT  =  KT  /  WT:  FV  =  KV  /  WT:  GOTO  710 

630  INPUT  "longitudinal  displacement  factor  in  cgs  units";  FL 
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640  INPUT  "transverse  displacement  factor  in  cgs  units";  FT 

650  INPUT  "vertical  displacement  factor  in  cgs  units";  FV 

660  KL  =  FL  *  WT:  KT  =  FT  *  WT:  KV  =  FV  *  WT 

670  AS  =  "a" 

680  INPUT  "generate  a  magnetic  data  file  (y/n)";  AS 

690  IF  AS  =  "Y"  OR  AS  =  "y"  THEN  GOTO  710 

700  IF  AS  =  "N"  OR  A$  =  "n"  THEN  GOTO  850  ELSE  GOTO  680 

710  INPUT  "magnetic  data  file  name";  M$ 

720  ON  ERROR  GOTO  730:  GOTO  740 

730  RESUME  710 

740  OPEN  "O",  #1,  MS 

750  WRITE  #1,  LAT,  LNG,  DEC,  DIP,  HE,  PLM,  PTM,  PVM,  WT,  KL,  KT,  KV,  FL,  FT,  FV 

760  CLOSE 

770  GOTO  850 

780  INPUT  "magnetic  data  file  name";  MS 

790  ON  ERROR  GOTO  800:  GOTO  810 

800  RESUME  780 

810  OPEN  "I",  #1,  MS 

820  INPUT  #1,  LAT,  LNG,  DEC,  DIP,  HE,  PLM,  PTM,  PVM,  WT,  KL,  KT,  KV,  FL,  FT,  FV 

830  CLOSE 

840  DECR  =  DEC  *  CON:  DIPR  =  DIP  *  CON 

850  A$  =  "a" 

860  INPUT  "processing  data  entry  by  f i le/keyboard  (f/k)";  A$ 

870  IF  AS  =  »F"  OR  A$  =  "f"  THEN  GOTO  1450 

880  IF  A$  =  "K"  OR  A$  =  "k"  THEN  GOTO  890  ELSE  GOTO  860 

890  A$  =  "a" 

900  INPUT  "input  sampling  period  (y/n)";  AS 

910  IF  AS  =  "Y"  OR  A$  =  "y"  THEN  GOTO  930 

920  IF  AS  =  "N"  OR  AS  =  "n"  THEN  GOTO  960  ELSE  GOTO  900 

930  INPUT  "sampling  period  in  seconds";  DT 

940  IF  DT  <=  0  THEN  PRINT  :  PRINT  "must  be  greater  than  zero":  PRINT  :  GOTO  930 

950  GOTO  1090 

960  AS  =  "a" 

970  INPUT  "input  maximum  magnetic  signal  frequency  (y/n)";  AS 

980  IF  AS  =  "Y"  OR  AS  =  "y"  THEN  GOTO  1000 

990  IF  AS  =  "N"  OR  AS  =  "n"  THEN  GOTO  1030  ELSE  GOTO  970 

1000  INPUT  "maximum  magnetic  signal  frequency  in  Hertz";  MAXF 

1010  IF  MAXF  <=  0  THEN  PRINT  :  PRINT  "must  be  greater  than  zero":  PRINT  :  GOTO  1000 

1020  GOTO  1080 

1030  INPUT  "minimum  target  slant  range  at  CPA  in  meters";  MINR0 

1040  IF  MINR0  <=  0  THEN  PRINT  :  PRINT  "must  be  greater  than  zero":  PRINT  :  GOTO  1030 

1050  INPUT  "maximum  magnetometer  relative  speed  in  knots";  MAXVMK 

1060  MAXVM  =  MAXVMK  *  KON:  MAXF  =  2  *  MAXVM  /  MINR0 

1070  IF  MAXVM  <=  0  THEN  PRINT  :  PRINT  "must  be  greater  than  zero":  PRINT  :  GOTO  1050 

1080  DT  =  1  /  (2  *  MAXF):  REM  low  pass  filter,  Nyquist  sampling  rate 

1090  AS  =  "a" 

1100  INPUT  "input  integration  time  (y/n)";  AS 

1110  IF  AS  =  "Y"  OR  AS  =  "y"  THEN  GOTO  1130 

1120  IF  AS  =  "N"  OR  AS  =  "n"  THEN  GOTO  1210  ELSE  GOTO  1100 

1130  INPUT  "integration  time  in  seconds";  IT 

1140  IF  IT  >=  DT  THEN  GOTO  1170 

1150  PRINT  :  PRINT  "IT  =  "  +  STRS(IT)  +  "  seconds     minimum  =  "  +  STRS(DT)  +  "  seconds":  PRINT 

1160  GOTO  1130 

1170  NS  =  2  *  INT< IT  /  DT  /  2)  +  1:  REM  adj  number  of  samples  per  integration  time 

1180  IF  NS  <=  151  THEN  GOTO  1340 

1190  PRINT  :  PRINT  "IT  =  "  +  STRS(IT)  +  "  seconds     maximum  =  "  +  STR$(150  *  DT)  +  "  seconds":  PRINT 

1200  GOTO  1130 

1210  INPUT  "maximum  target  slant  range  at  CPA  in  meters";  MAXR0 

1220  IF  MAXR0  <=  0  THEN  PRINT  :  PRINT  "must  be  greater  than  zero":  PRINT  :  GOTO  1210 

1230  INPUT  "minimum  magnetometer  relative  speed  in  knots";  MINVMK 

1240  MINVM  =  MINVMK  *  KON 

1250  IF  MINVM  <=  0  THEN  PRINT  :  PRINT  "must  be  greater  than  zero":  PRINT  :  GOTO  1230 

1260  IT  =  2  *  MAXR0  /  MINVM 

1270  IF  IT  >=  DT  THEN  GOTO  1300 

1280  PRINT  :  PRINT  "IT  =  »  +  STRS(IT)  +  "  seconds     minimum  =  "  +  STRS(DT)  +  "  seconds":  PRINT 
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1290  GOTO  1100 

1300  NS  =  2  *  INTCIT  /  DT  /  2)  +  1:  REM  adjusted  number  of  samples  per  integration  time 

1310  IF  NS  <=  151  THEN  GOTO  1340 

1320  PRINT  :  PRINT  "IT  =  "  +  STR$(IT)  +  "  seconds     maximum  =  "  +  STR$(150  *  DT)  +  "  seconds":  PRINT 

1330  GOTO  1100 

1340  A$  =  "a" 

1350  INPUT  "generate  a  processing  data  file  (y/n)";  A$ 

1360  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  1380 

1370  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  1510  ELSE  GOTO  1350 

1380  INPUT  "processing  data  file  name";  P$ 

1390  ON  ERROR  GOTO  1400:  GOTO  1410 

1400  RESUME  1380 

1410  OPEN  "O",  #1,  P* 

1420  WRITE  #1,  DT,  IT,  NS 

1430  CLOSE 

1440  GOTO  1510 

1450  INPUT  "processing  data  file  name";  P$ 

1460  ON  ERROR  GOTO  1470:  GOTO  1480 

1470  RESUME  1450 

1480  OPEN  "I",  #1,  P$ 

1490  INPUT  #1,  DT,  IT,  NS 

1500  CLOSE 

1510  A$  =  "a" 

1520  INPUT  "kinematic  data  entry  by  file/ keyboard  (f/k)";  A$ 

1530  IF  A$  =  "F"  OR  A$  =  "f"  THEN  GOTO  1790 

1540  IF  A*  =  "K"  OR  A$  =  "k"  THEN  GOTO  1550  ELSE  GOTO  1510 

1550  INPUT  "magnetometer  course  in  decimal  degrees  (0  if  at  rest)";  CM 

1560  INPUT  "magnetometer  speed  in  knots";  VMK 

1570  INPUT  "magnetometer  altitude  in  meters  (below  0  is  -)";  AM 

1580  INPUT  "target  course  in  decimal  degrees";  CT 

1590  INPUT  "target  speed  in  knots";  VTK 

1600  INPUT  "target  depth  in  meters  (above  0  is  -)";  AT 

1610  CMR  =  CM  *  CON:  CTR  =  CT  *  CON 

1620  WXK  =  VMK  *  SIN(CMR)  -  VTK  *  SIN(CTR) 

1630  WYK  =  VMK  *  COS(CMR)  -  VTK  *  COS(CTR) 

1640  Z  =  AM  +  AT:  REM  vertical  separation  (-  for  magnetometer  below  target) 

1650  X  =  WXK:  Y  =  WYK:  GOSUB  4170 

1660  CR  =  T:  C  =  CR  /  CON:  REM  magnetometer  relative  course 

1670  WK  =  R:  REM  magnetometer  relative  speed 

1680  A$  =  "a" 

1690  INPUT  "generate  a  kinematic  data  file  (y/n)";  A$ 

1700  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  1720 

1710  IF  AS  =  "N"  OR  A$  =  "n"  THEN  GOTO  1860  ELSE  GOTO  1690 

1720  INPUT  "kinematic  data  file  name";  K$ 

1730  ON  ERROR  GOTO  1740:  GOTO  1750 

1740  RESUME  1720 

1750  OPEN  "0",  #1,  K% 

1760  WRITE  #1,  CM,  VMK,  AM,  CT,  VTK,  AT,  Z,  C,  WK 

1770  CLOSE 

1780  GOTO  1860 

1790  INPUT  "kinematic  data  file  name";  K$ 

1800  ON  ERROR  GOTO  1810:  GOTO  1820 

1810  RESUME  1790 

1820  OPEN  "I",  #1,  K$ 

1830  INPUT  #1,  CM,  VMK,  AM,  CT,  VTK,  AT,  Z,  C,  WK 

1840  CLOSE 

1850  CMR  =  CM  *  CON:  CTR  =  CT  *  CON:  CR  =  C  *  CON 

1860  MCMR  =  (CMR  -  DECR):  REM  magnetometer  magnetic  course  in  radians 

1870  MCTR  =  (CTR  -  DECR):  REM  target  magnetic  course  in  radians 

1880  ILM  =  KL  *  HE  *  COS(DIPR)  *  COS(MCTR):  ITM  =  -KT  *  HE  *  COS(DIPR)  *  SIN(MCTR) 

1890  IVM  =  KV  *  HE  *  SIN(DIPR) 

1900  DMX  =  (PLM  +  ILM)  *  SIN(MCTR)  +  (PTM  +  ITM)  *  COS(MCTR) 

1910  DMY  =  (PLM  +  ILM)  *  COS(MCTR)  -  (PTM  +  ITM)  *  SIN(MCTR) 

1920  DMV  =  PVM  +  IVM 

1930  X  =  DMX:  Y  =  DMY 
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1940  GOSUB  4170 

1950  OMLR  =  T:  REM  dipole  azimuth  relative  to  magnetic  north 

1960  X  =  DMV:  Y  =  R 

1970  GOSUB  4170 

1980  DM  =  R:  OMR  =  T:  REM  dipole  depression  angle 

1990  OML  =  OMLR  /CON:  OM  =  OMR  /  CON 

2000  A$  =  "a" 

2010  INPUT  "generate  a  combined  magnetic,  processing  &  kinematic  data  file  (y/n)";  A$ 

2020  IF  AS  =  "Y"  OR  A$  =  "y"  THEN  GOTO  2040 

2030  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  2210  ELSE  GOTO  2010 

2040  INPUT  "combined  magnetic,  processing  &  kinematic  data  file  name";  E$ 

2050  ON  ERROR  GOTO  2060:  GOTO  2070 

2060  RESUME  2040 

2070  OPEN  "O",  #1,  E$ 

2080  WRITE  #1,  LAT,  LNG,  DEC,  DIP,  HE,  PLM,  PTM,  PVM,  WT,  KL,  KT,  KV,  FL,  FT,  FV,  DT,  IT 

2090  WRITE  #1,  NS,  CM,  VMK,  AM,  CT,  VTK,  AT,  Z,  C,  WK,  ILM,  ITM,  IVM,  DM,  OML,  OM,  M$,  P$,  K$ 

2100  CLOSE 

2110  GOTO  2210 

2120  INPUT  "combined  magnetic,  processing  &  kinematic  data  file  name";  E$ 

2130  ON  ERROR  GOTO  2140:  GOTO  2150 

2140  RESUME  2120 

2150  OPEN  "I",  #1,  E$ 

2160  INPUT  #1,  LAT,  LNG,  DEC,  DIP,  HE,  PLM,  PTM,  PVM,  WT,  KL,  KT,  KV,  FL,  FT,  FV,  DT,  IT 

2170  INPUT  #1,  NS,  CM,  VMK,  AM,  CT,  VTK,  AT,  Z,  C,  WK,  ILM,  ITM,  IVM,  DM,  OML,  OM,  M$,  P$,  K% 

2180  CLOSE 

2190  DECR  =  DEC  *  CON:  DIPR  =  DIP  *  CON 

2200  CR  =  C  *  CON:  OMLR  =  OML  *  CON:  OMR  =  OM  *  CON 

2210  INPUT  "required  false  alarm  rate  in  false  alarms  per  hour";  FAR 

2220  PF  =  FAR  *  IT  /  3600:  REM  false  alarm  probability 

2230  Y  =  PF:  IF  PF  >  .5  THEN  Y  =  1  -  Y:  REM  inverse  normal  approximation 

2240  Y  =  SQR(LOG(1  /  Y  /  Y>) 

2250  Y  =  Y  -  (11  +  Y  *  (12  +  13  *  Y))  /  (1  +  Y  *  (14  +  Y  *  (15  +  16  *  Y))) 

2260  IF  PF  <  .5  THEN  Y  =  -Y 

2270  ZP  =  -Y 

2280  CHI  =  NS  *  (1  -  2  /  9  /  NS  +  ZP  *  SQR(2  /  9  /  NS))  A  3:  REM  inverse  chi -square  approximation 

2290  INPUT  "magnetic  noise  in  gamma";  SIG 

2300  INPUT  "maximum  lateral  range  in  meters";  LRM 

2310  INPUT  "lateral  range  step  in  meters";  ST 

2320  IF  ST  <=  LRM  THEN  GOTO  2340 

2330  PRINT  :  PRINT  "maximum  step  is  "  +  STRS(LRM)  +  "  meters":  PRINT  :  GOTO  2310 

2340  NC  =  2  *  INT(LRM  /  ST)  +  1:  REM  number  of  lateral  range  function  values 

2350  IF  NC  <=  71  THEN  GOTO  2370 

2360  PRINT  :  PRINT  "minimum  step  is  "  +  STR$(LRM  /  35)  +  "  meters":  PRINT  :  GOTO  2310 

2370  AIT  =  DT  *  NS:  REM  adjusted  integration  time 

2380  W  =  WK  *  KON:  REM  magnetometer  relative  speed  in  meters/second 

2390  DS  =  W  *  DT:  REM  distance  between  samples  on  the  relative  track 

2400  X0  =  -(NC  -  1)  /  2  *  ST 

2410  FOR  I  =  0  TO  NC  -  1 

2420  X0(I)  =  X0 

2430  X  =  X0:  Y  =  Z 

2440  GOSUB  4170 

2450  R0  =  R:  R0(I)  =  R:  REM  target  slant  range  at  CPA  in  meters 

2460  DELR  =  T:  REM  target  depression  angle  complement  at  CPA  in  radians 

2470  IF  R0  =  0  THEN  GOTO  2800:  REM  zero  lateral  range  f»nd  vertical  separation  2480  DMF  =  DM  /  10  /  R0  *  3: 

REM  dipole  moment  factor 

2490  CMR  =  CR  -  DECR:  REM  target  relative  magnetic  course 

2500  L  =  COS(OMR)  *  COS(CMR  -  OMLR) 

2510  M  =  COS(DELR)  *  COS(OMR)  *  SIN(CMR  -  OMLR)  -  SIN(DELR)  *  SIN(OMR) 

2520  N  =  -SIN(DELR)  *  COS(OMR)  *  SIN(CMR  -  OMLR)  -  COS(DELR)  *  SIN(OMR) 

2530  L1  =  COS(DIPR)  *  COS(CMR) 

2540  M1  =  COS(DELR)  *  COS(DIPR)  *  SIN(CMR)  -  SIN(DELR)  *  SIN(DIPR) 

2550  N1  =  -SIN(DELR)  *  COS(DIPR)  *  SIN(CMR)  -  COS(DELR)  *  SIN(DIPR) 

2560  A2  =  2  *  L  *  L1  -  M  *  M1  -  N  *  N1:  REM  Anderson  Function  Coefficient 

2570  A1  =  3  *  (N  *  L1  +  L  *  N1):  REM  Anderson  Function  Coefficient 

2580  A0  =  2  *  N  *  N1  -  L  *  L1  -  M  *  M1:  REM  Anderson  Function  Coefficient 
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2590  SUM  =  0:  HMAX  =  0:  HMIN  =  0 

2600  FOR  J  =  0  TO  NS  -  1 

2610  S  =  (J  -  (NS  -  1)  /  2)  *  DS:  BA  =  S  /  RO:  REM  Anderson  Function  Argument  2620  AF  =  1  /  (1  ♦  BA  *  BA) 

2.5:  REM  Anderson  Function  Factor 

2630  HSF  =  (A2  *  BA  *  BA  +  A1  *  BA  +  A0)  *  AF:  REM  magnetic  signal  factor 

2640  HSU,  J)  =  DMF  *  HSF:  REM  magnetic  signal  value 

2650  IF  HS(I,  J)  >  HMAX  THEN  HMAX  =  HS(I,  J) 

2660  IF  HS(I,  J)  <  HMIN  THEN  HMIN  =  HS(I,  J) 

2670  SUM  =  SUM  +  HSF  *  HSF 

2680  NEXT  J 

2690  H(I)  =  HMAX  -  HMIN 

2700  K(I)  =  SQR(SUM) 

2710  W  =  -ZP  +  DMF  *  SQR(SUM)  /  SIG 

2720  LAM  =  DMF  *  DMF  *  SUM  /  (SIG  *  SIG):  AA  =  NS  +  LAM:  BB  =  1  +  LAM  /  (NS  +  LAM) 

2730  ZN  =  -SQR(2  *  CHI  /  BB)  +  SQR(2  *  AA  /  BB  -  1):  X1  =  W 

2740  GOSUB  4230 

2750  IF  Y1  >  1  THEN  Y1  =  1 

2760  PDCC(I)  =  Y1:  X1  =  ZN 

2770  GOSUB  4230 

2780  IF  Y1  >  1  THEN  Y1  =  1 

2790  PDSL(I)  =  Y1 

2800  XO  =  XO  +  ST 

2810  NEXT  I 

2820  C  =  C  /  360:  C  =  (C  -  INT(C))  *  360 

2830  IF  C  <  0  THEN  C  =  360  +  C 

2840  OML  =  (OML  +  DEC)  /  360:  OML  =  (OML  -  INT(OML))  *  360 

2850  IF  OML  <  0  THEN  OML  =  360  +  OML 

2860  PRINT  :  A$  =  "a" 

2870  INPUT  "generate  a  program  data  file  (y/n)";  A$ 

2880  IF  AS  =  "Y"  OR  A$  =  "y"  THEN  GOTO  2900 

2890  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  3230  ELSE  GOTO  2870 

2900  INPUT  "program  data  file  name";  D$ 

2910  ON  ERROR  GOTO  2920:  GOTO  2930 

2920  RESUME  2900 

2930  OPEN  "0",  #1,  D$ 

2940  WRITE  #1,  LAT,  LNG,  DEC,  DIP,  HE,  PLM,  PTM,  PVM,  WT,  KL,  KT,  KV,  FL,  FT,  FV,  DT,  IT,  AIT 

2950  WRITE  #1,  NS,  CM,  VMK,  AM,  CT,  VTK,  AT,  ILM,  ITM,  IVM,  Z,  C,  WK,  DM,  OML,  OM,  FAR,  PF,  SIG,  ST 

2960  WRITE  #1,  LRM,  DS,  NC,  M$,  P$,  K$,  E$ 

2970  FOR  I  =  0  TO  NC  -  1 

2980  WRITE  #1,  X0(I),  PDCC(I),  PDSL(I),  K(I),  H(I),  R0(I) 

2990  NEXT  I 

3000  FOR  I  =  0  TO  NC  -  1 

3010  FOR  J  =  0  TO  NS  -  1 

3020  WRITE  #1,  HS(I,  J) 

3030  NEXT  J 

3040  NEXT  I 

3050  CLOSE 

3060  GOTO  3230 

3070  INPUT  "program  data  file  name";  D$ 

3080  ON  ERROR  GOTO  3090:  GOTO  3100 

3090  RESUME  3070 

3100  OPEN  "I",  #1,  D$ 

3110  INPUT  #1,  LAT,  LNG,  DEC,  DIP,  HE,  PLM,  PTM,  PVM,  WT,  KL,  KT,  KV,  FL,  FT,  FV,  DT,  IT,  AIT 

3120  INPUT  #1,  NS,  CM,  VMK,  AM,  CT,  VTK,  AT,  ILM,  ITM,  IVM,  Z,  C,  WK,  DM,  OML,  OM,  FAR,  PF,  SIG,  ST 

3130  INPUT  #1,  LRM,  DS,  NC,  M$,  P$,  K$,  E$ 

3140  FOR  I  =  0  TO  NC  -  1 

3150  INPUT  #1,  X0(I),  PDCC(I),  PDSL(I),  K(I),  H(I),  R0(I) 

3160  NEXT  I 

3170  FOR  I  =  0  TO  NC  -  1 

3180  FOR  J  =  0  TO  NS  -  1 

3190  INPUT  #1,  HS(I,  J) 

3200  NEXT  J 

3210  NEXT  I 

3220  CLOSE 

3230  PRINT  :  A$  =  "a" 

58 


3240 

INPUT 

3250 

IF  A$ 

3260 

IF  A$ 

3270 

LPRINT 

3280 

LPRINT 

3290 

LPRINT 

3300 

LPRINT 

3310 

LPRINT 

3320 

LPRINT 

3330 

LPRINT 

3340 

LPRINT 

3350 

LPRINT 

3360 

LPRINT 

3370 

LPRINT 

3380 

LPRINT 

3390 

LPRINT 

3400 

LPRINT 

3410 

LPRINT 

3420 

LPRINT 

3430 

LPRINT 

3440 

LPRINT 

3450 

LPRINT 

3460 

LPRINT 

3470 

LPRINT 

3480 

LPRINT 

3490 

LPRINT 

3500 

LPRINT 

3510 

LPRINT 

3520 

LPRINT 

3530 

LPRINT 

3540 

LPRINT 

3550 

LPRINT 

3560 

LPRINT 

3570 

LPRINT 

3580 

LPRINT 

3590 

LPRINT 

3600 

LPRINT 

3610 

LPRINT 

3620 

LPRINT 

3630 

LPRINT 

3640 

LPRINT 

3650 

LPRINT 

3660 

LPRINT 

3670 

LPRINT 

3680 

LPRINT 

3690 

LPRINT 

3700 

LPRINT 

3710 

LPRINT 

3720 

LPRINT 

3730 

LPRINT 

3740 

LPRINT 

3750 

FOR  I 

3760 

LPRINT 

3770 

NEXT  I 

3780 

PRINT 

3790 

INPUT 

3800 

IF  A$ 

3810 

IF  A$ 

3820 

LPRINT 

3830 

LPRINT 

3840 

LPRINT 

3850 

LPRINT 

3860 

LPRINT 

3870 

FOR  I 

3880 

LPRINT 

"print  encounter  parameter  values  (y/n)";  AS 

=  "Y"  OR  AS  =  "y"  THEN  GOTO  3270 

=  "N"  OR  A$  =  "n"  THEN  GOTO  3780  ELSE  GOTO  3240 

"program  file  name 
"program  data  file  name 
"magnetic  data  file  name 
"processing  data  file  name 
"kinematic  data  file  name 

"combined  magnetic,  processing  &  kinematic  data  file  name 

"encounter  latitude  (decimal  degrees)  " 

"encounter  longitude  (decimal  degrees)  " 

"encounter  variation  (decimal  degrees)  " 

"encounter  dip  angle  (decimal  degrees)  " 

"encounter  magnetic  field  intensity  (gamma)  " 

"permanent  longitudinal  moment  (oersted-cm3)  " 

"permanent  transverse  moment  (oersted-cm3)  " 

"permanent  vertical  moment  (oersted-cm3)  " 

"target  displacement  (tons)  " 

"target  longitudinal  permeabi I ity  coefficient  " 

"target  transverse  permeability  coefficient  " 

"target  vertical  permeability  coefficient  " 

"target  longitudinal  permeability  factor  " 

"target  transverse  permeability  factor  " 

"target  vertical  permeability  factor  " 

"sampling  period  (seconds)  " 

"integration  time  (seconds)           '  " 

"adjusted  integration  time  (seconds)  " 

"number  of  samples  per  encounter  " 

"magnetometer  course  (decimal  degrees)  " 

"magnetometer  speed  (knots)  " 

"magnetometer  altitude  (meters)  " 

"target  course  (decimal  degrees)  " 

"target  speed  (knots)  " 

"target  depth  (meters)  " 

"magnetometer  relative  course  (decimal  degrees)  " 

"magnetometer  relative  speed  (knots)  " 

"magnetometer- target  vertical  separation  (meters)  " 
"target  induced  longitudinal  dipole  moment  (oersted-cm3)" 
"target  induced  transverse  dipole  moment  (oersted-cm3)  " 
"target  induced  vertical  dipole  moment  (oersted-cm3)    " 

"magnetic  dipole  moment  (oersted-cm3)  " 

"dipole  moment  azimuth  (decimal  degrees)  " 

"dipole  moment  depression  angle  (decimal  degrees)  " 
"distance  between  samples  on  the  relative  track  (meters)" 

"false  alarm  rate  (false  alarms  per  hour)  " 

"false  alarm  probability  " 

"magnetic  noise  (gamma)  " 

"maximum  lateral  range  (meters)  " 

"lateral  range  step  (meters)  " 

"number  of  lateral  range  function  values  " 
=  0  TO  15 


"  + 


"  + 


N$ 
D$ 
M$ 
P$ 
K$ 
1  +  E$ 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 


:  A$  =  "a" 

"print  lateral  range  function  values  (y/n)";  A$ 

=  "Y"  OR  A$  =  "y"  THEN  GOTO  3820 

=  "N"  OR  A$  =  "n"  THEN  GOTO  3900  ELSE  GOTO  3790 

D$;  "  lateral  range  function  values" 

:  LPRINT 

"L         p(cc)       p(sl)        H  R 

"meters  gamma        meters" 


LAT 

LNG 

DEC 

DIP 

HE 

PLM 

PTM 

PVM 

WT 

KL 

KT 

tcv 

FL 

FT 

FV 

DT 

IT 

AIT 

NS 

CM 

VMK 

AM 

CT 

VTK 

AT 

C 

UK 

Z 

ILM 

ITM 

IVM 

DM 

OML 

OM 

DS 

FAR 

PF 

SIG 

LRM 

ST 

NC 


k" 


=  0  TO  NC  -  1 
X0(I);  TABOO);  PDCC(I);  TAB(24);  PDSL(I);  TAB(38);  H(I);  TAB(52);  R0(I);  TAB(70);  K(I) 
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3890  NEXT  I 

3900  PRINT  :  A$  =  "a" 

3910  INPUT  "print  magnetic  signal  values  (y/n)";  A$ 

3920  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  3940 

3930  IF  AS  =  "N"  OR  AS  =  "n"  THEN  GOTO  4120  ELSE  GOTO  3910 

3940  PRINT  :  PRINT  "identify  signal  by  encounter  lateral  range  index" 

3950  PRINT  "lateral  range  equals  the  index  times  "  +  STRS(ST)  +  "  meters" 

3960  PRINT  "index  values:  -"  +  STR$((NC  -  1)  /  2)  +  "  to  "  +  STR$((NC  -  1)  /  2) 

3970  INPUT  "lateral  range  index";  K:  I  =  K  +  (NC  -  1)  /  2 

3980  LPRINT  :  LPRINT  :  LPRINT 

3990  LPRINT  D$;  "  magnetic  signal  values" 

4000  LPRINT  "for  a  lateral  range  of  ";  K  *  ST;  "  meters" 

4010  LPRINT 

4020  LPRINT  "relative  CPA  distance";  TAB(35);  "magnetic  signal" 

4030  LPRINT  "s  in  meters";  TAB(35);  "Hs  in  gamma" 

4040  LPRINT 

4050  FOR  J  =  0  TO  NS  -  1 

4060  LPRINT  (J  -  (NS  -  1 )  /  2)  *  DS;  TAB(35);  HS(I,  J) 

4070  NEXT  J 

4080  PRINT  :  A$  =  "a" 

4090  INPUT  "print  magnetic  signal  values  for  a  different  lateral  range  (y/n)";  AS 

4100  IF  AS  =  "Y"  OR  AS  =  "y"  THEN  GOTO  3940 

4110  IF  AS  =  "N"  OR  A$  =  "n"  THEN  GOTO  4120  ELSE  GOTO  4090 

4120  PRINT  :  A$  =  "a" 

4130  INPUT  "continue  to  use  the  program  (y/n)";  A$ 

4140  IF  AS  =  "Y"  OR  A$  =  "y"  THEN  GOTO  10 

4150  IF  AS  =  "N"  OR  AS  =  »n"  THEN  GOTO  4160  ELSE  GOTO  4130 

4160  END 

4170  R  =  SQR(X  *  X  +  Y  *  Y):  REM  rectangular  to  polar  conversion 

4180  IF  R  =  0  THEN  T  =  0:  RETURN 

4190  IF  ABS(X  /  R)  =  1  THEN  Q  =  SGN(X)  *  (PI  /  2)  ELSE  Q  =  ATN(X  /  R  /  SQR(1  -  X  *  X  /  R  /  R)) 

4200  IF  ABS(Y  /  R)  =  1  THEN  T  =  (PI  /  2)  *  (1  -  SGN(Y))  ELSE  T  =  (PI  /  2)  -  ATN(Y  /  R  /  SQR(1  -  Y  *  Y  /  R  / 

R)) 

4210  IF  Q  <  0  THEN  T  =  2  *  PI  -  T 

4220  RETURN 

4230  Y1  =  X1:  IF  X1  <  0  THEN  Y1  =  -Y1:  REM  normal  approximation 

4240  G  =  1  /  (1  +  Q0  *  YD 

4250  Y1  =  EXP(-Y1  *  Y1  /  2)  /  SOR(2  *  PI )  *  G  *  (01  +  G  *  (02  +  G  *  (03  +  G  *  (04  +  G  *  05)))) 

4260  IF  X1  >=  0  THEN  Y1  =  1  -  Y1 

4270  RETURN 
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